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The flow of a compressible fluid with weak entropy 


changes 


By ROY GUNDERSEN 
Department of Mathematics, Illinois Institute of Technology 


(Received 7 August 1957 and in revised form 27 November 1957) 


SUMMARY 

Perturbations of a given flow are considered, and the equations 
which govern the one-dimensional, non-steady flow of an inviscid, 
ideal compressible gas are linearized in the neighbourhood of this 
known solution, assumed isentropic, by a formal perturbation 
expansion. ‘The ‘perturbed flow is not assumed isentropic. 
Explicit solutions are obtained for a basic flow which is uniform 
or a centred simple wave and for an arbitrary simple wave if 
the perturbed flow is isentropic. 

The perturbation of a uniform shock and perturbations in a 
shock tube lead to functional equations of a particular type, and 
a discussion of their solution is given. 

A similar analysis is used to discuss the flow in a tube of slowly 
varying cross-section. 


1. INTRODUCTION 

The one-dimensional, non-steady flow of an inviscid compressible gas 
which, furthermore, is assumed to be ideal and polytropic with constant 
specific heats is governed by a system of non-linear, first-order partial 
differential equations which is always hyperbolic. Given an isentropic 
flow, a solution of the governing equations, this paper will discuss 
‘perturbations’ of this basic flow; i.e. this flow is disturbed by some 
means, and it is desired to find the resultant flow. 

To linearize in the neighbourhood of this known solution, the dependent 
variables are formally expanded as functions of a small parameter 6, viz. 

t) = u,(x, t) + du,(x, t) + d7uo(x, t)+..., 
, t) = C(x, t) + dc, (x, t) + 67c,(x, t) +..., 
, t) = 59+ 65,(x, t) + 6°s,(x, t) + 

where u, c and s are respectively the particle velocity, the local speed 
sound and the specific entropy, and the zero subscript denotes the basic 
(known) flow. When these expansions are substituted into the governing 
equations and grouped with respect to increasing powers of 6, the terms 
free of 6 are satisfied identically, so that the resulting equations may be 
divided through by 6 and, taking the limit as 6 > 0, the equations governing 
the first approximation to the behaviour of the disturbed flow are obtained. 
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‘These equations are linear with the same characteristic surfaces as the 
original system and form a system of non-homogeneous first-order partial 
differential equations. It is convenient to attack the problem by first 
considering the associated homogeneous system, which has the physical 
significance that, to the order of terms retained (i.e. to the order of 5), the 
perturbed flow ts still isentropic. ‘The problem is consequently reduced 
to the determination of a particular solution of the non-homogeneous 
equations. Since the equations for the perturbations are linear, super- 
position of flows is possible. 

Explicit solutions are obtained for a basic flow which is uniform or a 
centred simple wave, and for an arbitrary simple wave if the perturbed 
flow is isentropic. ‘The application of these results to the perturbation of 
a uniform shock and perturbations in a shock tube lead to functional 
equations of a particular type, and a discussion of their solution is given in 
Appendices I and II. The flow in a long tube whose cross-section is 
made to undergo a small variation with wx is treated by a similar analysis 
in part IIT. 

Many authors have discussed the application of perturbation techniques 
to gas-dynamical problems, but, in general, have linearized in the neighbour- 
hood of a constant state, with the result that the characteristics are rectilinear 
and parallel. For many flows, these characteristics are inadequate, and a 
more exact treatment of the characteristics is necessary. In the present 
paper, the question of correcting characteristics is not met, for the exact 
characteristics of the initial flow are always employed. A summary of 
$2, §3 and $5 has already been given by Germain & Gundersen (1955). 


PART I. GENERAL THEORY 


2. "THE BASIC EQUATIONS AND SOLUTION FOR AN ORIGINALLY UNIFORM FLOW 


The equation of state is 


S—s™ 
P=exp| —— |p’; 
Cy 


and the equations which characterize the one-dimensional non-steady flow 


are 


where P is the pressure, p the density, s* the specific entropy at some 
reference state, y the ratio of the specific heats at constant pressure c, 
and at constant volume c,, and the subscripts denote partial derivatives. 


v 


Throughout, it is assumed that an isentropic solution is known, namely, 
Uy(x, t), Co(x, t), Sy = constant. A formal linearization in the neighbourhood 


\ 
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of this known solution gives the following equations which govern the first 
approximation to the behaviour of the disturbed flow: 


Cyt ily Ce t My — Leg ty, +, C9, + Hy — 1)tto, €, = 9, (2.1) 


Ky —1)uy + My —1)tty ty + C0 Cr + My — 1 )ttg Hy £1 Cor = CRS 7/2Cpy, — ( 


/ 


Sap+Ug Sy, = 9. (Z.3) 


These equations are linear with the same characteristic surfaces as the 
original system. 

According to (2.3), s; remains constant along the particle paths of the 
given flow, i.e. along dx/dt = uy. Also, po(dx —u, dt) is the exact differential 
of a function 4%, which, when equated to a constant, defines the particle 
paths. Consequently, the solution of (2.3) is 


$1 = wy), (2.4) 


with w an arbitrary function. It is convenient to define a new function 


H,(x, t) by 


é Siz = ee w’' (ho) = vy = l)c, Hy. (2.5) 

‘Therefore, the two equations which serve to determine u,, c, are 
Cut Uo Cre t3(y — 1 )cg yy + Uy Co, + Hy — 1)g, Cc, = 9, (2.6) 
Uy + Ug Uy. + 2(y — 1)~eg Cy, +g Uy + 2(y — 1)~'e9,.¢, = Ay. (2.7) 


It is convenient to introduce the characteristic parameters of the basic 
How 
ply + Co/(y—1) = «, — 3Ug + C/(y—1) = B, 


and the functions 


pty +¢,/(y—1) = A, —}u,+¢,/(y—1) = B. 
The combinations (2.6)/(y — 1) + (2.7) give 
A, + (uy +¢,)A,+ S[A(y +1) +(y —3)B] x, = SA, (2.8) 
B, + (ug — &)B,+ HAGB-y)—(v + 1)BIB, = —tAh. (2.9) 


\ perturbation in which the perturbed flow is still isentropic corresponds 
to a solution of the homogeneous (H, = 0) system associated with the 
linear system (2.8) and (2.9). 

For the case of an initially uniform flow, uw) and cy are constants, and 


therefore x, = 8, = 0. ‘The homogeneous system then admits the following 
general solution, in terms of two arbitrary functions F and G of one argument 


A = Fl[x—(uy+cy)t], B = G[x — (uy —cy)t]. (2.10) 


In the non-isentropic case, from (2.4) and (2.5), Hy is a function of 
~—u,t, and 

Hy = 2x'(x— wot), CoA = 6p B= x(x—tgt), (2.11) 

where y is an arbitrary function, is a particular solution for which u, = 0; 

i.e. the addition of an entropy perturbation affects c,(x, t) but not u,(x, ft). 


2N2 
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3. PERTURBATION OF A SIMPLE WAVE : GENERAL THEORY 

For the case of the simple wave, one of the characteristic parameters 
of the basic flow remains constant, e.g. B = By. B is then determined 
from (2.9) and A from (2.8). ‘The solution for the perturbation of an 
arbitrary simple wave will be given but, first, the centred simple wave 
will be discussed in some detail. 

The wave is assumed to be centred at the origin and characterized by 
> { 


8 = By. From the definition of the characteristic parameters, 


Uy = %~— Bo, Chis 3(y 1)(« Li Bo); 


? , me 
MoE ge 
y+1{¢t : ed | 


Substituting these results into (2.8) and (2.9) specialized to the case 
of an isentropic perturbed flow, the following equations are obtained: 


tA,+xA,4 | 4 a B | = (3.1) 
4 


“4 a oy 1 
th, T E 7 1 orm Bo ‘|B, = Q). (3.2) 
4 V 4 


B will be determined first by solving (3.2), and, after this result has been 
substituted into (3.1), 4 may easily be determined. ‘The characteristics 


therefore 





R 
| 
| 








of (3.2) may be written as 
dt yt tds (3 + ) ol ae (3.3) 
t (3 —y)x—4(y —1) Bot x+2Pyt ’ 
and dB =(). It follows from (3.3) that 
t= K,(x+2B,t)’tve™ 


with A, a constant; and since 


oy 1 e+ 2 Bot 
Co = are 1 t aod ’ 


the representation of the curvilinear characteristics of the simple wave 1s 


t K3(Cy t)” 1)/(3—y) 








3 
eT, 
: / 








with A, a constant, or 
pytat = const. (3.4) 
This form is not classical, but will be seen to be a very convenient one. 
If F is an arbitrary differentiable function, it 1s convenient to write the 


solution of (3.2) as 


B = (py co)!2t F'(py Co #2). (3.5) 
Consequently, the solution of (3.1) is 
3-y 1 . 
At — = F (py Cy t?) + G(x/t), (3.6) 


yl pay 


where G is an arbitrary function. 
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The particle paths are given by pg cyt = const., or 


y= = const. 
For the non-isentropic case, it is convenient to write 
LH, = (2c, t)w,(y). 
The following particular integral (to be added to the previous solution) is 
easily obtained : 


uy 0, cy = iy t 1)(¢9/v)o,(y), 
with w, an arbitrary function and, as in the perturbation of a constant state, 
the addition of the entropy perturbation affects c,(x,f) but not w,(x, ¢). 
This question will be taken up in detail after the arbitrary simple wave 
has been discussed. 
For the arbitrary simple wave (isentropic perturbed flow), let x (2), 
t,(2) be the parametric representation of a curvilinear characteristic of the 


given flow. ‘lhe wave may be represented by 
x Yo(2) u [uo(2) u Co(3) 7, t ty(3) rks (3.7) 
5Ug(2) Co(2)/( 1) 30 (3.8) 


Equation (2.9) reduces to 
B,+(u g—cy)B, = 9, (3.9) 
i.e. B is constant along any curvilinear characteristic. ‘lhe characteristics 
of (3.9) are dB = 0 and dx = (u,—c,) dt, or, from (3.7), 
dXy d dr dt, 


po —_ 9 . “ 
dz rs d= (Uy T Co) T dz aCy (Uy 6) a ’ 
hence, as (X»),¢,)) is a curvilinear characteristic, the following differential 
) ) 
equation is obtained for 7(z): 


2c us Eu : (Uy +C)) = 0 
= "dz dz sd 0 
or, utilizing (3.3), 
l dr ly+11 de, 
rdz 29-16 dz” 


This has the solution 


2 


PolyT? = const., OF poCy(¢t—t,)* = const. 

This result for the representation of the curvilinear characteristics of 
the simple wave could have been written down immediately by analogy 
with (3.4). ‘The solution for B is 

b 2(p9 Co)! *(t ty) F’| py eo(t ty)”] (3.10) 


in terms of an arbitrary, differentiable function F. ‘The characteristics 


of (2.8) may be written as 
dx 2dA 


ute [3-y)B-(y+1)A]o,’ 


dt (3.11) 
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lhe integral of the first equation of (3.11) has the first family of character- 
istics (rectilinear characteristics) for level curves. Let this first integral be 
3(x, ¢) = const. (3.12) 

This function is given implicitly by 

X— Xo(%) = [uo(z) + Co(2)][t — to(2)]- 
From the definition of the characteristic parameters, it follows that 

ko 2(y t 1) [Uo rly ay 3)]Bo- 


Substitution in (3.11) yields the following differential equation for A: 


2¢,(dty dz)y 11dA 3 Y > 
E — ty on “(dey/dz) | Tt 4 A ary RB. 


To integrate this equation, B is substituted from (3.10) and the first integral 
(3.12) is employed, i.e. s may be considered as constant. ‘This yields 
another first integral, and the solution for A(t, 2) is 


2¢)(dt,/dz)y—1 
aie ee” be oe UA 
[ fo (dc,/dz) y+ |4 2) 
3-y 1 : 


at a Fool —49)") (3.13) 
Y T Po ©o és 


in terms of the arbitrary functions F and G. In the non-isentropic case, 





the problem is reduced to two ordinary differential equations, but it does 
not seem possible to solve these two equations explicitly in the general 
case. ‘The question of actually determining these arbitrary functions in a 
specific problem will be discussed later 


4. CRITERION FOR THE EXISTENCE OF A PARTICULAR SOLUTION w= 0 IN THE 
NON-ISENTROPIC CASE 
In the non-isentropic perturbation of a constant state and of a centred 
simple wave, the addition of the entropy perturbation just affects the sound 
speed, whereas the particle velocity is the same for isentropic and non- 
isentropic perturbed flows, i.e. there exists a particular solution u, = 0. 
Unfortunately, such a nice result does not obtain for the arbitrary simple 
wave. It is of interest to determine under what conditions such a result 
is to be expected, for the answer is not at all obvious. 
It is equivalent to look for a solution of (2.6) and (2.7) with u, = 0, viz. 
Cyp + Uy Cy, + 3(y — 1g, Cy = 9, (4.1) 
(C1 Cox = (Po 5/2, yw (Yo), (4.2) 
dibs, = po(dx — Up dt). 
Thus, it is clear that the problem consists of expressing the compatibility 
of equations (4.1) and (4.2), and the corresponding identities will determine 
which functions u(x, t), ¢o(x, t) allow this compatibility. 
The characteristics of (4.1) are 
dx 2dc, 
Uy i. (y-1 Won 4 


dt = 
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Irom the first two ratios, a first integral is 
dx — Uy dat 0, or Wj = const. 
A second one will be found easily since co(x,t) is a particular solution 
of (4.1). This first integral is c,/cg = const. Hence, the solution of (4.1) is 
C9 iF (ho) (4.3) 
where f is an arbitrary function. Substitution of (4.3) into (4.2) gives 
Cos | w' (uo) 
Po Lo af 2c, 
Hence, the necessary and sufficient conditions for the existence of a 
particular solution wu, = 0 is 


Cor = ] 


( Cy, = a function of zp. 


Po ©o 
By taking the material derivative, this condition becomes 
) (). 

Consequently, in order to have a particular solution with u, = 0, it ts 
necessary and sufficient that uo(x, t)is a linear function of x. Then cy = cy f(b) 
with a convenient function f. \t is necessary to investigate all solutions of 
this form, i.e. vu) = a(t)x+(t). ‘This has been carried out but will not be 


included here. 


5. AN APPROXIMATION FOR FLOWS WHICH DIFFER ONLY SLIGHTLY FROM 
\ UNIFORM STATE 

In general, it does not seem possible to solve (2.8) and (2.9) explicitly 
for an arbitrary initial flow. However, the preceding techniques may be 
utilized if the given flow does not differ greatly from a uniform state or a 
simple wave. 

For example, consider the problem of a flow, which is nearly a uniform 
flow (to, Cg). Denote by u,, ¢, s; and ty, Cy, So the terms of first and second 
order obtained by a formal perturbation method. Let 

My = Ug tty, E=Cytey tls GyH Cot ey, B= UgtUyt+ug, S$ = $,+5p. 
By termwise addition of the equations determining u,, ¢, 5; and Ug, Cy, Sy 
and neglecting terms of third order (e.g. Cy, Uy C,), the following equations 
are obtained: 

C+ ty €,+ Hy —1)e, a, = 9, 
Ui, + ty, UW, +2(y--1)1G,€, = E75,./v(y (5.1) 
5, + ti, §, = 0. 
Formally, this system of equations is identical to (2.1), (2.2) and (2.3) if 
the terms u,,. and cy, are neglected. Consequently, the following principle 
has been established. Jf i, and ¢, define the flow to the second order, ti and é 
define the flow to the third order. 

The solution of the system (5.1) is immediate. With the introduction 

of the quantities A, B, x and f£, the system may be written as 


Lain &. 1 Hy tys 
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and may be solved by quadratures. In the isentropic case, A = F(a), 
B = G(8) where F and G are two arbitrary functions. ‘hus, it is clear 
that the analysis is exactly that which is met in the classical theory (the 
Euler— Poisson equation). 

Note that the neglect of the terms up», and cg, does not affect the 
principal part of each differential equation, so that the characteristics 


remain unchanged. 


6. ISENTROPIC PERTURBATION OF AN ARBITRARY SIMPLE WAVE: 
DETERMINATION OF ARBITRARY FUNCTIONS 
In $3, the perturbation of a simple wave was solved in terms of two 
arbitrary functions. ‘The present and subsequent section are devoted to the 
question of the determination of these arbitrary functions. It will be shown 
that, at least theoretically, no difficulties are met. 


t 


x=X(t) | 


—— A - 
x *. aa 
eis 
\ as ow x=CRt ae 
> a 
* ie Sa 
ho, nx 


Figure 1. 


Consider a tube of gas, initially at rest with constant sound speed c,, 
and constant density p,, from which a piston, originally at rest, is withdrawn 
with increasing speed until a constant final velocity is attained. A simple 
wave, which advances with the speed c,, into the gas at rest, is generated. 
‘he (x, ¢)-representation of the flow is shown in figure 1. The region 
c,t, denoted in figure 1 by the symbol R, is a constant state with 
(u,c) = (0,c,). The flow due to the receding motion of the piston is 
restricted to the region x < c,t, and the region denoted by O is a simple 
wave. As the piston attains its final velocity, the flow becomes a constant 
state, denoted by the symbol 2, with particle velocity equal to the piston 
This is a classical problem, and the solution is well known. 


X 


velocity. 
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It is now assumed that a small perturbation takes place and that the 
resultant flow is still isentropic. ‘his problem was solved in §3 in terms 
of two arbitrary functions. ‘The discussion of this section will be in rather 
general terms without concern for computational details, for it is the basic 
principles, and not the details, which are important here; they will provide 
the basis for the solution of the perturbation of the centred simple wave. 

Let the following initial conditions be given for t = 0: 

u = 02(x), ( Cp t+dh(x), 
where g and / are known functions. 
From the perturbation of a constant state (§2), the solution for x > c,t 


is 
in = n(x—Cpt)—E(x+c,p b), 2(y --1) 1é, = o(x—ept)+&(xt+e, t), 
where 
n(x) = be(x)+h(x)/(y—1), E(x) 58(x) + A(x)/(y — 1), 
and perturbations are denoted by a bar. Also, let the piston velocity be 


u = % = X(t)+ SX(2). 


‘lhe rectilinear characteristics of the simple wave are given by = = const., 
and the curvilinear cross-characteristics by 


9 


Po Colt — to)” = 4, 
where a@ is a constant. Recall that [x9(=), fo(s)] is the parametric repre- 
sentation of a non-linear characteristic of the given simple wave. 

‘lhe function z was defined implicitly, and it is assumed that this relation 
has been inverted to give s(x,t). ‘This value is substituted into p,(z), etc., 
which gives functions which will be denoted by pp,(x,¢), etc. ‘These are 
substituted into (3.10) and (3.13), and the resultant functions are called 
B = B,(x,t), A = A(x, t), G(z) = G,(x,2) and F,,. 

The perturbations introduced at ¢ = 0 travel along the characteristics 
in region FR to the first characteristic of the simple wave, viz. x = c,t. 
At a point (x,,¢,) = (¢p¢,,¢,) of this characteristic, the perturbations are 

A,(x,,t,) = 7(9), B, (x1, t)) = &(2cp ty). (6.1) 


This perturbation (a%;,¢,;) travels along the non-linear characteristic 
through the point (x,,¢,;) to the piston curve. Let the equation of this 


non-linear characteristic be 


Por Colt — ty)” = 1, (6.2) 
and let the intersection of (6.2) with the piston path, v = X(t), be (xj, ¢)). 
The conditions (6.1) on « = c¢,t allow the determination of /, from 

2( P91 Cor )"'2(t — to) F [01 Cor (¢ — tor)*] = €(2€R t1)- (6.3) 


Thus F,, and therefore B,. is known throughout the simple wave region. 
In general, 4 = 4,—B,; therefore, at the point (x;,¢;|) of the piston 
curve, we have 


X(t) = A,(x!, t}) — B,(, #1). (6.4) 
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vy — 1 (dt dz) ] 
a 9 Y—! \Mo/4@2) , 
D! [ ty “OU +1 (decd) | in 


As the first and last terms in (6.4) are known, this gives a condition for the 
determination of G,(x}, tj): 
G,(xt, ti) = Di[X,(t}) + B,(x}, t1)] 4 
ot F [pt c(t —211)?]. (6.5) 
y +1 (phi cin)” = sas 7 
Clearly, this analysis applies at each point of the first characteristic of the 
simple wave, i.e. at each point the perturbation travels along the non-linear 
characteristic through the point to the piston curve. Consequently, (6.5) 
holds, with appropriate changes in the arguments, all along the portion 
of the piston curve in the simple wave region and allows the determination 
of G,(x,t). From classical theory, this is a well-set problem (in Picard’s 
terminology, a third problem) as the value of B, is known along the 
characteristic v = c,t and the value of A, is known along the piston curve 


Now let 


which has time-like orientation. 

‘The function G, determines the perturbations which travel along the 
rectilinear characteristics of the simple wave. ‘lhe perturbed flow in 
region 2 is determined by data along the last characteristic of the simple 
wave and on the rectilinear portion of the piston curve. 


7. ISENTROPIC PERTURBATION OF A CENTRED SIMPLE WAVE 
Consider the same original situation as in $6, but suppose now that the 
acceleration of the piston from rest to a constant terminal velocity wu, = —V 
(I’ > 0) takes place instantaneously. ‘lhe rectilinear characteristics of the 
rarefaction wave, generated by this piston motion, all pass through the 


origin as depicted in figure 2. 


Figure 2. 


On this known initial flow, the following initial conditions for the 
perturbed flow are superposed: t# =0, t =0; &, = (y—I1)f(x), v > 0; 
u(0,0) =0; ¢(0,0) = (y—1)f(0). The analysis of $6 will be used to 
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determine the solution of the present problem. ‘The subsequent analysis 
will also give information about the value of é,(0, 0). 

The utilization of the solution for the perturbation of a constant state 
and initial conditions gives 


A=f(x-cet) B=f(x+ext). (7.1) 


The application of (7.1) on x = c, 


G(cr) = f(0)t + —-——~5 F (pp cpt? 
(Cp) F(0) y4 i (ppc, yr lr ent ); 
f(2cpt) = App cp)*t F'(pp ext). 


In the limit as ¢t > 0, 


t and substitution in (3.5) and (3.6) gives 





one ] 
G > = ——_—_ *(( 
(a) = 1 peg | (7.2) 
f(9) = (), 


where it is assumed that f is continuous and that F’ is bounded near ¢ = 0. 
This last statement implies that the perturbations are bounded near ¢ = 0. 

Equation (7.2) says that, for the present problem, it is necessary that 
€(0,0) = 0, i.e. both perturbation velocities must vanish at the origin 
where the basic flow is discontinuous. 

The present problem is clearly a limiting case of the problem of § 6. 
‘The G-function, a function of x/t, determines the perturbations which 
travel along the rectilinear characteristics of the simple wave. It is 
determined by the perturbations carried to the piston curve along the 
curvilinear cross-characteristics of the simple wave. In the limit as the 
rectilinear characteristics form a pencil of lines through the origin, the 
G-function will be determined by the perturbations at the origin only. 
As the perturbations are zero at the origin, it is clear that G is at most a 
constant, and this constant will be chosen as zero for convenience. 

It is easily found that 


F(p3¢3t?) = (prep) | “f(2cn¥) dy. 


The solution in the simple wave is now known, and the perturbed flow 
in region | is determined by conditions along the last characteristic of the 
simple wave and along the piston. 


PART II. APPLICATIONS 
As a first application, the perturbation of a shock produced by a uniform 
compressive motion of a piston will be considered. A piston, initially at 
rest, is at ¢ = 0 suddenly pushed with constant velocity into a gas which 
is also initially at rest. A uniform shock is produced, which moves with 
constant velocity into the gas, and the gas behind the shock is in a steady 
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state with gas velocity equal to the piston velocity. ‘The state of rest in 
front of the shock is denoted by R, and the constant state behind the shock 
is denoted by 0. ‘The shock velocity is given by 
Wo = (y+ 1)uy + [eh + {Cy + 1)ty/4}7]"?. 

Suppose that the piston is given a small prescribed perturbation for ¢ > fy. 
This will perturb the shock locus, and, as this is the perturbation of a strong 
shock, the flow behind the shock will no longer be isentropic; thus the 
entropy perturbations must be considered. It is convenient to solve this 
problem in two steps. First, the inverse problem of the determination 
of the perturbed piston path when the shock is given a small perturbation 
is solved, and, second, this solution is used to solve the direct problem. 


8. ‘THE INVERSE PROBLEM 

Suppose that, from some particular time fj, the shock is given a small 
prescribed perturbation; i.e. let the shock velocity be W = W, for t <t, 
and WW = W,+<(t) for t > t., where «(t) is a known function sufficiently 
small so that terms of higher order than the first may be neglected. ‘This 
perturbation will cause a corresponding perturbation of the particle velocity ; 
let the new particle velocity near the shock be u = uy+7(t), where x(t) is 
also small and must be a function of ¢(t) as will be shown. 

‘The perturbations of tp, ¢y and sy, which will be denoted by a, ¢) and Sy, 
will be found in terms of ¢(t) along the shock. With these boundary 
conditions, the solution for the non-isentropic perturbation of a constant 
state may be utilized to find the perturbed flow behind the shock. ‘Then, 
fiy(uyt, t) gives the perturbation of the piston velocity. 

It is convenient to tabulate many of the standard relations between 
the flow parameters in perturbation form. Appendix III contains some 
of the more important relations. 

By the use of the Prandtl relation in perturbation form (Appendix ITI, 
relation VIII), it follows that 

iy = 2K, hh), Cy = (y—1)K, €(t), 
where 


Ky = $[2(1-0—up/Wy], Ky = 2p +1) (Wy — up) [ey + 03/2 Wo ey, 


Consequently, 
A =(K,+K,)e(t) = Ky e(é), B =(K,—K,)e(t)  K,e(t) (8.1) 


on the shock. 
From relation V of Appendix III, 


Py u,—(1—@)W, 
Po Keo ey i 


and from relation I of Appendix ITI, 
(2c, , (y 1) ; [up - (1 6)W,]) - ) 
0 < Co (y 1 )Ky VW Cy ” _— re(t) K. e(t). (8.2) 


BY 
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In general, from (2.11) and (2.5), 
Sy = 2cyv(y— leg “x(x — uy t) 


which implies, from (8.2), that 


Ke “xX —Upt oe 
X(% — Uy t) 5 Pe DS 7... (8.3) 


throughout region O. 

By the use of the general solution for the perturbation of a constant 
state and the conditions (8.1) and (8.3), and applying the boundary conditions 
on the unperturbed shock locus*, the following solution is obtained: 


y cy K; X—Ugt , Cy K; xX — (Uy + Co)t 
cies ae yy - 5‘ a : = | | K, 2c, v(y - 5 | a7 Uy — Al 
_— C) Kz xX—Uyt ; by Far xX — (Uy — Cp)t 
sine 2c¢, (y ia a Wy, | | K, CoVly "es ~ Ug t ef 


be ' ta K, V— (Up + Co)t 
il [ | 2eyyx(y- ake - ee 
Cy K; V — (Uy — Cy)t 
"S| ee ee a al a §. 
| es 1) Ky |¢| ar Uy 4 | a) 


is the particle velocity perturbation throughout region O, and, in particular 
its value on the piston path is obtained by setting v = ugt in (8.4). The 
perturbed piston and shock paths are obtained by simple integration. 

















bo 








Figure 3. 





It is not immediately obvious how the characteristics should be oriented. 
However, in figure 3, consider an arbitrary point on the piston curve 


*'The error thus introduced is of the order of the neglected terms and therefore 
neglible, viz. if the original shock locus is x W, t, the new shock locus is 


x Wot (' e(y) dy, and 


= . *} 
n| Wo t- |, €(V) dy| (Wo t) +7’( Ws t) | e(y) dy 
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(e.g. ?*), and let t_ be the intersection of the shock locus and the character- 
istic with slope u)—cy through this point. The characteristic with slope 
ly + Cy through ¢* intersects the shock locus in a point which will be denoted 
by ¢.. ‘These points are given by 

t Co t* |(ug + Co — Wo), t_ = Cyt*/(cy—uy t+ W,). 
Consequently, from (8.4) with v = uot, it is seen that a perturbation at a 
particular point of the piston curve is determined by the perturbation at 
two points of the shock curve with the depicted orientation. Hence, the 
(x, t)-representation is as in figure 4. The region in which the flow is 
perturbed is cross-hatched. ‘The relation between f, and f, is 


ty = Cyt,/(ug t+ co— W,).- 


aN " 
] pas 
A 
| Ly 
L 
k 
k i > al 
Le > 
rf t> “i 
A YJ 
Le fe P 
i> 
<b 
L ux n+ 
| t, w De = 
x 
é = 
Figure 4. 


9. DIRECT PROBLEM 
From the solution of the inverse problem, the solution of the direct 
problem is easily obtained, and the solution is given by a uniformly 
convergent infinite series. Let 7(t) be the prescribed perturbation of the 
piston velocity, and ¢(t) the perturbation of the shock velocity. From (8.4), 


we get 


y(t) K.—--a 


[ Cy K; y Col ] 
| K = 
| 2c, rN) L) F W, Uy + Cy 


| 
= 


f ‘| rf T 
“o ° | “0 1 
Rye caesar Reel ei 4. £91) 
1 ; 2 | 


| Co — Uy + Wy | Uy +Cy— W, 


Figure 5 illustrates the scheme to be utilized. ‘The lines drawn between 
the piston and shock paths are the characteristics, the odd-numbered lines 


having slope (uw) +c¢)) and the even-numbered lines having slope (tu — cy). 
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Now, 7 evaluated at the point 1 is given in terms of multiples of € evaluated 
at the points 1’ and 2’. (2) is given in terms of multiples of ¢(2’) and «(3’), 
etc. «(1’) will be determined. Symbolically, take the first-mentioned 
relation (1, 1’, 2’), equation (9.2), and for 2’ substitute its equivalent from 
equation (9.3). ‘This gives a relation (1, 3, 1’, 4’). By continuing this 
process, a relation of the form (1, 3, 5, 7, ..., ”, 1’, [u+1]’) for nm odd is 
obtained. By indefinite continuation of this process, a series representation 
is obtained for € in terms of multiples of 7 evaluated at the points K, where 
AK = 1, 3,5,.... ‘The general term of this series is determined and proved 
by induction, and the convergence is established. 


| . 9 
Po. 2 
/ & 

| / 3 pe 

| Te / 

| / Se 3! 

| 3K6 rd 

1g D4 


aan 3 


Figure 5. 





The coordinates of the point 1’ are chosen to be (x, ¢) = (W,?*, ¢*), 
and the following notation will be employed: 


Uy tey-- Wy = 7, W—uyteo = 4 Cot* = t,. 
By induction, 
y(1) = R, €(rt,/g?) + Ry €(t,/r) (9.2) 
n(2) = R, €(r?t,/¢°) + Ry €(t,/q) 
(3) = R, €(r3t,/q*) + Re €(rt,/q") (9.3) 
n(n) = R, €(r"t,/q"*1) + Ry €(r"*t,/q"—). (9.4) 


From (9.2), 
e(t,/r) = n(1)/R.—(R,/Ro)e(rt,/ 9"). 
‘The last term is replaced by its equivalent from (9.3), with the result 
e(t,/r) = n(1)/R.—(R,/R5)y(3) + (R,/Re)*e(77t,/9"). 


For the first m terms, the series is 
e(t,/r) = > (- 1)A RK RSF y(2K — 1) +(— 1)"(R,/Ro)"e(r?"“1t, /g@").. (9.5) 


Mathematical induction is used to prove this assertion, which is certainly 
true for K = 1. Assume it is true for K =n. From (9.4), it follows that 
fey, g”) y(2n +1) R,—(R, R,)e( 2 tlt /q™ 2). (9.6) 
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Substitution of (9.6) into (9.5) gives 
I 


e(t,/r) = > o(2K—1)(—1)4 RE UR, B+ (—1)"4(R,/ Ry)! the(r2” +12, 1g?" +2), 
K=1 
(9.7) 
‘Therefore, the assumption that (9.5) is true for 7 implies it is true for n+ 1. 
Hence, (9.5) is true for all n. 

Since it is assumed that all perturbations are small, the assumption is 
made that « and y are uniformly bounded. With this reasonable assumption 
and the easily-seen fact that R,/R, < 1, it is readily established that, as 
n> oo, the series is uniformly convergent. Hence 

e(t,/r) = = f 1)4 ‘Re 'R, *y(2K—1). 
K-1 
Functional equations of the type (9.1) are treated from a general point of 


view in Appendix I. 


10. PERTURBATIONS IN A SHOCK TUBE 

The orthodox shock tube consists of a tube divided by a membrane 
(say at x =Q0) which separates two fluids, not necessarily the same, 
originally at rest but not at the same pressure. ‘his is one frequently used 
device for generating constant shocks in a long tube. It will be assumed 
that the fluid in the region «x 0 is at the higher pressure. If, at ¢ = 0, 
the membrane is removed instantaneously, the fluid in the region x > 0 
will expand, while the fluid in the region x < 0 will be compressed. Under 
these assumptions, the expansion takes place through a centred rarefaction 
wave, and the compression across a uniform shock. In general, there will 
be a contact discontinuity which separates the two portions of fluid as 


depicted in figure 6. 


At 











Figure 6, 


On this basic flow, which is assumed to be completely known, a 


perturbation is imposed. ‘The perturbation conditions are not assumed 
to be the same for the two portions of the fluid. In particular, the isentropic 
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character of the initial flow will not be retained for the region x > 0. Let 
these conditions be: for t=0 and x> 0, u =dg(x), c=c¢,+6A(x), 
s=5)+5f(x); for t=0 and x <0, u=5n(x), c=c.+dm(x), s = 5. 
Here c,, c. and sy are constants, and g, h, f, m and n are known functions. 
For continuity at x = 0, it is assumed that g(0) = n(0), A(O) = m(0) and 
f(0) = 0. ‘The contact discontinuity is given by x = Kt. 

It is not possible to solve independently for the perturbed flows on the 
two sides of the contact discontinuity. For the region x > 0, no difficulty 
is encountered in regions 1 and 3. The solution in 3 gives conditions on 
the last characteristic of the simple wave, but these are insufficient to 
determine the flow in region 5. However, if on the contact discontinuity 
i; is specified as €,(t), which is not known, the flow in region 5 can be 
easily determined in terms of &,(t). For the region x < 0, the flow will 
be non-isentropic since the shock is perturbed. By expressing a, ¢, and 5, 
in terms of the shock velocity perturbation W, the flow in region 4 can be 
determined in terms of W, which is not known. 

Across the contact discontinuity, there is continuity of pressure and 
particle velocity, and it is easy to see that this implies continuity of the 
corresponding perturbations. ‘This continuity condition yields two 
equations in €,(t) and W(t). By elimination of €,(t), a functional equation, 
treated in detail in Appendix II, is obtained for W and, once W has been 
determined, the solution of the problem is easily completed. 

In terms of (x) and &(x), as defined in §6, the solution in region 1 is 


ty = (x—C,t)— E(x +e, t)— 3p. f'(x- Cyt) + dpi f(xt+e; t), 
2(y —1)-1, = n(x—c, t) + E(x +c, t)—3p,¢, f'(x—c,t)- 
— 3p, 0, f(x +c, t) +p, ¢ f(x), 


ry 


5 = vv — lever f(*). 

These perturbations travel along the characteristics in region 1 to the 
first characteristic of the centred simple wave, where they supply boundary 
conditions for determining the perturbed flow in region 3. On x= ¢,f, 
the first characteristic of the simple wave is 

— . 1 ane ti . 1 Sey (7 
B = &(2c, t)— 3p, ¢, f’(2c, t) + dp, ¢, f'(ce, 8), 
Fin = Wy — Vey pr f(r *)- 

Application of these conditions and the solution for the non-isentropic 

perturbation of a centred simple wave gives 


2 8 (eslexy!? 9067 
w4(egt") = Ac p, | yf ey) dy, 
“(Ps CalPrer)’ ¢ : 
F(ps¢3t?)/(p,&)'? = | 1§(2¢, 3) — 3p, Cy f'(2¢, 2) + 3p a f'(4.3)} de— 
—p,ci| 2? | v'f'(e,y) dydz. 
J 0 y=0 


F.M. 20 
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Consequently, 4, B and §, are known throughout the simple wave. 
As the expressions are very complicated, let the values on the last 
characteristic of the simple wave be 

A= C,(¢), B = C(t), 55, = C3(¢), 


where ¢,, ¢ and @, are known functions. However, these conditions are 
not sufficient to determine the perturbed flow in region 5. ‘This is clear, 
for it cannot be expected that data along one characteristic, but nothing 
on another curve, would yield a well-set problem. It is convenient to let 
u; = €,(t) on x = Kt, the contact discontinuity, where, of course, €,(t) ts 
not known at the present time. ‘The perturbed flow in region 5 can now 
be easily found. Simple reasoning shows that ¢, is a constant multiple 
of €,(t) plus and a known function and, therefore, 


P, = vé(t)+ T(t) 


where v is a constant and |" is a known function. 

As the pressure and velocity are continuous across the contact dis- 
continuity, the perturbations of these quantities are also continuous across 
this surface, so that 

i= 6(t) Py =vé (t+ 
on x = At. 

For region 2, the solution is easily found. ‘This gives # and ¢, and, 
from Appendix IT], all perturbed quantities in region 2 are easily determined. 

From Appendix III, the values of the perturbed quantities in region 4 
may be found just behind the shock in terms of W. The analysis is simple 
though tedious, and simple reasoning shows that they have the following 
form: 

P, = R;, +R, W, p, = R, +R, W, é,= R, +R, W,) 
Sy Ry, i Ri» W, My Ris ' Ri Vv, ‘ (10.1) 
A=R,;,+R,, W, B= Ry + Ris W, 


where all the R’s are known, and odd-numbered subscripts denote functions 
of t while even-numbered subscripts denote constants. 

Utilization of the solution for the non-isentropic perturbation of a 
constant state and the boundary conditions (10.1) gives 


X=Ugt " (u,+c,)t 
Rul | , Ral Woe | 


| x—(uUy—C,)t z,[ ¥—uyt 
-R., Hl ro | R,, H lw | 


where all the R’s are known. 
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Throughout the region 4, we have 


P, . 2yp, Z P, : 
aly—1) *” Ge, 
The continuity of P across the contact discontinuity gives an equation ¢ 
the form 


K—-u,- 
E(t) = Ro;(t) + Roo W | (7 a) + 
ty 4 
| K—-uy+e, U4 
aL are) |e [Cra 00 


and the continuity of @ across the contact discontinuity gives 


é a "7 tug — Cy ,| (K Uy tq 
E,(t) = Rop(t) + Reg W | or =) |; Ry W | (7 — ad (10.. 


Division of (10.2) by (10.3) serves to eliminate €,(t), with the result 


B K-—uyt+e, K —4,—¢, 
Ryo W| (aera) |+R, W ash 


where all the R’s are known. Once W has been found, €,(t) is given 





o>) 
— 


immediately, and the solution is easily completed. ‘This functional equation 
is discussed in Appendix LI. 


PART III. TUBE OF SLOWLY VARYING CROSS-SECTION 

‘The flow in a long tube whose cross-section is made to undergo a small 
variation with x can be treated by a similar analysis. It is convenient to 
express the cross-section of the tube in the form E(x) = E, +6 F(x) 
where FE, is the original uniform cross-sectional area and 6 is a small 
parameter. 

It will be assumed that the tube consists of two parts, one of uniform 
cross-section, and one of perturbed cross-section. A uniform hana 
introduced in the tube where the cross-section is uniform, will travel along 
the tube and remain uniform until it reaches the transition section. ‘lhe 
gas in front of the shock is assumed to be at rest, and the small effects due 
to the discontinuity in the slope of the tube at the transition section will 
be neglected. After the shock passes the transition section, it wili be 
perturbed, and the flow behind the shock will necessarily be non-isentropic. 
The problem is to describe the flow behind the shock and to find the 
perturbed shock locus. The results are applied to several problems, and 
some comparisons are made with work by a previous author. 


11. DERIVATION OF EQUATIONS WHICH GOVERN THE FLOW 
Clearly, the Euler and energy equations are not attected by the cross- 
sectional perturbation, but the continuity equation contains an additional 
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term. ‘lhe mass, contained in an arbitrary volume 7, may be written as 
M | pdr = | p E(x) dx. 


The continuity equation, DMW/Dt = 0, imphies that 
(pE), + (pEu), = 0. 
Consequently, the basic equations are 
u,+uu,+2c,/(y—1) =s,c/e,y(y—1), 
2Ec,|(y — 1) + 2Euc,|(y —1)+cEu,+ck,+uck, = 0, 
5, +us, = V. 
A formal linearization in the neighbourhood of an originally isentropic 
flow gives 
Up Uy, + Uy + 2Cy Cel (y — 1) + Uy Up + 2egy y/(y — 1) = C517 /€py(y — 1), (11-1) 
Qe uly — 1) + €p Wy, + 29 Cy, /(y — 1) + €1 Up, + 2u, Co,-/(y — 1) 4 
+ Uy Cy Ey ./Eg + 2¢y E,/Ey(y — 1) + €p U9, Fi /Eo 4 
+ 2tty Coy Ey/EgQ(y—1) = 0, (11.2) 
Sy + Up $1, = 9. (11.3) 


Y 


It is thus clear that the cross-sectional perturbation merely introduces 
a non-homogeneous term into the basic perturbation equations which were 
previously derived, so that the mathematical results of the previous sections 
may be utilized, and the problem is reduced to the determination of a 
particular integral of (11.1), (11.2) and (11.3). 

For the case of an initially uniform flow, (11.3) may be solved 
independently, and, utilizing the same procedure as employed in §2, 
the following equations are obtained: 

A,+(uUpt%)A, - bug Cy Ey ,/Eo 4 
B, + (ty — C9) B, = — $9 Co E1,/Eo — 5 Ho. 


‘H, 


The general solution in terms of three arbitrary functions is 
H, = 2y'(x— uot), 
A = Flx—(uy + Cot] + ¢y 1x(% — Uo t) — $tlg Co Ey /Eg(g +o), (11.4) 
B = G[x — (ug — Cp)t] + Cy x(x — Ug t) — Sg Cy Ey /Eg(to— co). (11.5) 
Consequently, there is a perturbation, due to the entropy variations, 
which travels along the particle paths, and this is measured by y; the 
perturbations measured by F travel along one family of characteristics 
with sonic velocity relative to the fluid, and the perturbations measured 
by G travel along the other family of characteristics with sonic velocity 
relative to the fluid. 
Chester (1953) has considered the disturbance produced behind a 
plane shock of arbitrary strength travelling down a two-dimensional channei 
of non-uniform width, and the problem was linearized on the basis of small 


variations of the width of the channel. ‘lhe variations in width were assumed 
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to take place within a finite length of the channel, and this region of transition 
separates two uniform portions each of infinite length but not necessarily 
equal widths. ‘lhe pressure field behind the shock was built up from the 
known solution of the diffraction of a shock wave travelling along a wall 
with a corner (Lighthill 1948). In a later paper, Chester (1954) considered 
essentially the same problem, but a different approach was used and applied 
to a tube of arbitrary cross-section. 

Let the flow in front of the shock be given by m = 0, cy, Po, po, the 
How behind by wy, cz, Ps, pz, and the shock velocity by w. Let v. = u,—w 
and M, = w/e). ‘The following relations are easily obtained: 


u, = 2w(1— M, “)/(y +1), 
Ve U,—-w = wWy—1+2M,*)/(y +1), 


M, = \Vo9\/¢o = [(y —1+2M, *)/(2y 


/ 


m = U/C, = 2(1— MP)/M,(y + 1). 


~ 


1)M,, *)]#2, 


Chester restricted his treatment to the investigation of the average 
pressure and, for the case where the tube consists of two cylinders connected 
by a transition section, the following result was obtained: 

Py =KMP,- Py )E,/Eo, 

K* = 2(1+m)(1+M,7+2M,)". 
‘(he parameter A* decreases monotonically with the shock strength, and 
0:5 > K* > 0-394. Actually, it is clear that the disturbance depends 
only on the variations in the area of the tube and not on the actual shape 
of the cross-section, as will be seen shortly. 

It will now be shown that, using the purely one-dimensional analysis 
of this section, Chester’s result can be obtained quite simply. Consider 
the general solution (11.4) and (11.5). ‘The term in F£, is due directly to 
the changes in cross-section, and this disturbance is reflected at the shock 
and gives rise to the term involving G. From the way Chester has 
formulated the problem (i.e. with the shock coming from infinity, so to 
speak), the contribution given by F is not included, for there is no mechanism 
(e.g. a piston curve) for reflection upstream of the shock which could give 
rise to such aterm. In other words, Chester’s result should be obtainable 
from (11.4) and (11.5) by putting # = 0. ‘The details are sketched below. 
If F = 0, it follows from (11.4) that 
m(1 +m) 'E, 


2c, ¢(y—1) 2¢,y(y—-1) ZE, 
Krom the shock perturbation tables given in Appendix III, it is easily 
determined that 
Ga|ealy — 1) — $igleyr(y — 1) = 4P aly Py 
Hence 


thy /Cy P,/yP,—mE,/(1+m)E,. (11.6) 


Zia 
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Us/Cy = 2(M, — M,, 1)/(y +1), (y + 1)iig/cg = 201+ MZ 2)M,, 
and, from relation LX of Appendix II], 
P,/P, = 4yMy My/(2yMi.—y +1), 








then uw (14 M?)(2yM2—y+1)P, 
Co 2y(y + 1)M3 P, ; 
cs Pspy — (2yMi—yv + Il2+(v—-1)Mi) 
i tia || 
Hence c; Pops Mey + 1) 


Cy ly|Cy Cy = tay[Cy = (1 + M2)P,/2yM? M, P. 
Also, from (11.6), 
P1+M?+2M?M,)/2P, M2 = —yM,mE,/(1+m)E,. 


Finally, since 
P, = {2yM?—(y—1)}(P2— P,)/2y(M2- 1), 


then 


: (mM, (2yM2—y +1)) 
P, — —(P,—P,)B, E> 1-+m)-(1-+2M, + Mp2 {Maley MG y +P) 


Mj — 1 i. 
[t is easily shown that the term in braces is equal to 2. Hence, the result 
is the same as that obtained by Chester. 
‘lhe analysis of the present paper is useful because it can be applied 
to a variety of problems. In the next section, a problem where the 
perturbations are reflected on a piston curve is treated. ‘lhe analysis could 





be applied to a problem where the tube is open at one end. In such a 
problem, reflection at a surface of constant pressure would be encountered. 


12. SLOWLY CONVERGING OR DIVERGING CROSS-SECTION 
For this problem, E,(x) = Kx, where K > 0 for a diverging cross-section 
and K <0 for a converging cross-section. By a uniform compressive 
motion, a uniform shock is introduced at x = 0. It is assumed that the 


rf 


Figure 7. 
piston is pushed continuously with velocity V, and the original shock 
velocity is w. Figure 7 is the (x, t)-representation of the flow. Region O 
is a rest state, and the perturbed flow behind the shock is denoted by 2. 
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‘lhe problem will be solved by expressing the perturbations @,, ¢, and §, 
in terms of %, at present unknown, and using the general solution (§ 11) 
to find 4 and B throughout region 2 in terms of #. From the assumptions 
of the problem, # = 0 on the piston path; this gives a functional equation 
from which @ can be easily determined. 

The Prandtl relation, in perturbation form, gives 

i, = 2K, ®, é, = K,(y —1)w 

on the shock, with 


2K, = 2(1-4)-u 





; 2(w — Uy) us 2 m 
fs a EE oh es a = . : 
is Co(y +1) 2we, yr! 7 y+1 (] I, *) 


Consequently 
A =(K,+K,)# = K,2, B= (K,—K,)w = K,w 
on the shock. 
By the same analysis as in part II, 











26 KRaly (4 l)c 9 | A)w 5 
i 4 2¢, Kaly I) = iw i « : [Ms a ( eo] + w(t) 
Ce w(u., — w) 
¢(y—1) 4 (J Vf?) , y 1 —" 
~@ lpr ty +! oso = NY 


on the shock. ‘UVhroughout region 2 
§, = 2c, y(y — 1)ce? x(x — ust), 
so that 
2erpyV(y — 1)x(x— ust) = Ky co @[(x — uy t)/(w — uy)]. 
By the use of the general solution and the boundary conditions, the 


solution for region 2 is 








Uy Cy Kx -[| x—Ust 
A - ~—-———_ +. ¢, K.[2c,, y(y — 1)]} te] ——— |-4 
2E, (us + C9) 2 Ks ey W — Us 
f 2K. -[ x—(uy+Co)t 
< ee 2 Ks Se} — ~ (ue + Ca) | 
| "Qe, p(y — 1) WW — ts — €5 
; [x — (uy + €2)t] 
oe Kw en Sot... at. 
ae” [2E (us + 5) |[ — Us — Cs] 
Uy Cy Kx ()K -[ x—Uyt 
B a — — Ww ———— i 
2E,(u.—Cs) 2¢,y(y—1) Lw-us, 


xX—(Uy—Cs)t] usc, Ku 
W—Us+Cy | 2ZE(ug—C3) 


As uw =0 on x= Vt, the following equation is obtained for the 


determination of @: 


cK = Col = Col | rm e fa 
— | 0) — ——————— — 
| K, 2¢,.y(y » | i Us | E Us + =| y—1) i 
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It is easily seen that a linear function of ¢ satisfies this equation, and the 


solution is 


2c, v(y—1) 
‘The pressure perturbation behind the shock is, in terms of the parameters 
M,, M, and m, 


Eee ~K,—K, M, [28 (0) = Uy K(w — uy)t. 


P. mr K,(P2- Py )E, Ey , 


(y #1)K, = My 142M, | 20% +mM,(1—M, 2)+(1+M,°) 





(y ak? 2M; ) f 4(1 iS M?) Lt — 14 | I 
+m? —1+ ———~——, > | . 
2y 1G +1)(1—-M,?) (y—1+2M ?)} 
Here K, is a monotonically decreasing function of the shock strength; 
and, for y = 7/5, 0-259 < K, < 0-608. A graph of this function is given 


~) 


\ 


in figure 8. 
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Figure 8. 


13. CONCLUSIONS 
In this work, several problems of flows with weak entropy changes 
have been considered and solved explicitly. Without going further into 
the details of the problems treated, we may note the type of functional 
equations to which the solutions of these problems lead. ‘That is, the 
values of the function to be found appeared in the same relation for different 


values of the independent variable. 
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It is clear that in cases where an analytic solution cannot be found, the 
general equations are suitable for a numerical calculation developed step 
by step along the characteristics. ‘This procedure would be necessary for 
the perturbation of a flow with a shock of constant intensity. 

It is also clear that the problems treated here are among the most simple 
non-isentropic flows that can be considered. Undoubtedly, an analogous 
theory could be developed for steady flow in two dimensions. However, 
it is not at ail sure that analytic results as complete as the present ones 
could be obtained. Perturbations can also be considered which introduce 
a new coordinate, but the problem is then even more complex. Several 
particular cases where the perturbed flow is uniform have already been 
considered, but much important work remains to be done. 

The following problem, which has already been solved (Gundersen 
1954), may be noted. By the use of the method of characteristic coordinates, 
an analytic solution was obtained for the interaction, in an inviscid ideal 
gas with a constant ratio of specific heats of 5/3, of a centred rarefaction 
wave and a non-uniform shock of constant intensity, i.e. the entropy jump 
is constant all along the shock with the result that the flow is isentropic on 
both sides of the shock. ‘This work shows that the analysis of this paper 
need not be restricted to the perturbation of uniform shocks. 

These studies are important because flows with shock waves play an 
increasingly important role in numerous technical problems, and, except 
for certain exceptional cases, such flows are not isentropic. 


The problem discussed in this paper was suggested by Professor Paul 
Germain and the research directed by him. ‘lhe author wishes to express 
his sincere appreciation to Professor Germain for his valuable advice and 
his continuous and stimulating interest in the project. 


APPENDIX | 
‘The functional equation 


A f(at)+ Bf(bt) = g(t) (A 1) 
will be considered, where A, B, a and 6 are constants. Assume 
|B) : |Al, \b\ > Al, (A 2) 


and put 
alb=£, A/B=K, r=O5t, f(bt)= F(t), g(7/b) = BG(r), 
fat) = f(bta/b) = F (Er). 
Consequently, (A 1) may be written as 
K F(ér) + F(r) = G(r), f= hk, (A) <4. (A4) 
Conversely, if G(z) is known, and if F(z) satisfies (A 4), i.e. if (A 4) can be 
solved, (A 1) can also be solved. 


Theorem: If G(r) is continuous in an interval J (—N =< 7 < N), 
which implies that G(r) is also bounded in J, and if F(z) is uniformly bounded 
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in a neighbourhood of 7 = 0, there exists one and only one solution ot (A 4) 
defined in /. F(z) is a continuous function in /. If G(r) is continuously 
differentiable in J, F(r) is also continuously differentiable in J. 


Proof: lf vel, "7 El because & 1. ‘Then 
F(r) + K F(Er) = G(r), 
F(ér) + K F(&r) = G(ér), 


Consequently, 
F(r)+(—1)"K"*1F(E" +47) = G(r) + ai —1)"K?G(E€?r) = g,,(7). 
p=1 
Hence 
F(r)—g,(7)| = |K"*"|| F(E"+#7)]. 
As F(r) is to be uniformly bounded in a neighbourhood of 7 = 0, it is 


possible. to find n large enough so that F(€"*17) is bounded. Call this 
bound C. ‘Thus 


VN eee n+l Aur 
F(r) —g,(7)| < C|K"*"|. (A5) 
Since K 1, if such a solution F(r7) exists, it must be the uniform limit 
of the sequence g,(7) in J (uniform convergence). Conversely, if g,,(7) 


is a uniform sequence of functions in /, the limit is the unique solution 
of (A 4). However, it is obvious that the series 

> (—1)?K” G(&7) 

p ‘) 
is uniformly convergent in J. Hence the solution of (A 4) 1s 


F(r) = ¥ (—1)P APB? g(aPr/b? +), | 
p=0 | 


Remarks: 1. If G(r) = 7", F(z) T(1+ KE"), n 0. ‘Thus, if G(r) 


n 
may be approximated by } d,7?, F(z) may be approximated by 
p=1 
> 4, 7?/(1+ Ke). 
p=1 
[f G(r) has a ‘Taylor expansion in / (i.e. it is an analytic-function), F(z) is 
an analytic function in J. 
2. If K = 1, it is necessary to assume G(0) = 0. ‘To preserve unique- 
ness, it is necessary to assume f(t) is continuous near ¢ Q). 
3. If |é 1, there is no problem if €=1, K~—1; 
F(r) = G(r)/(1+ K). 
4. If é el 4 1, then 
K F(—7)+ F(z) = G(r), K F(r)+ F(—7) = G(—7), 
F(r)(1— K?) = G(r)+ K G(—7). 


He = 1,18 1, the problem is either impossible or indeterminate. 
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APPENDIX II 
‘The functional equation 


f(t) +hf(mt) + Kf(nt) = g(t) (B 1) 
will be considered, where h, A, m and nu are constants, and jm iy 


n| <1, [kA|+|K| < 1. 


Theorem: Uf f(t) is uniformly bounded in a neighbourhood of t = 0, 
and g(t) is continuous in an interval / (¢ A), there exists one and 
only one solution of (B1) defined in J. f(t) is a continuous function 
in / and, if g(t) is continuously differentiable in J, f(t) is also continuously 
differentiable in /. ‘he solution of (B1) is 


f(t) = > > (—1)? 4h? K4 g(mPnt)(p + gq)! /p!q!, (B 2) 
p=0q=0 
where this series is uniformly convergent. 
First, formally consider the functions g(m?n%t), where p,q = 0, 1, 2,3,..., 
and attempt a solution of (B1) of the form 


> Dd pq glmPn't). (B 3) 
Og= 
The series (B3) is a solution of (B1) if 
Coe l, €y4 7 Key o 15 C107 hey » QO, | 
, | (B4) 
Cn.a he Ly ow 1 0. J 
Put €y.g = (— 1) 4d, , WR“. 
‘Then 
Res d, Lg ae 1 doy = 1, dy diy I, 
and therefore 
d,, =(p+4q)!/p!¢q!. (B5) 


It is now necessary to show that (B2) is actually the solution of (B 1) 
and that it is unique. Consider the operator L, which is defined by 
Lp =X DY (— 1) hr K4 o(ment)(p +9) '/p'q), p+q<r. (B6) 
p= U @= a) 
The application of this operator to (B 1) gives 
, ee - (p+q) ! 
f(t)-L,gt) => D> APK hf(m?* nt) + Kf(mPntt)|(— 1)e*4 _ 
p+q=rvl pig: 
(2 +49)! 
es . 


! 


f-L.g) = > > LA f(met nat) +|K)\ f(mPnt*3t)\]|h|PK 
prq=rtl p:4q 
Since f is to be uniformly bounded in a neighbourhood of ¢ = 0, it is 
possible to tind r large enough so that 
f(mPnt+t)| < M, f(m?*'n4t)| < M, 
where M is a constant. ‘Therefore 
f-L.g M)> > hP\ KK “(pt+q)!/p!q!. (B 7) 
pt+q=r+l 


As r —- ow, the right-hand side approaches zero and 
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f( = La(t) y y (—1)? **h? K% 9(men’t)(p+q)!/p!q!. (BS) 
i) q ia] 


As g is uniformly bounded in /, say g(t) N, |t A, the series 


N[1+\Al+{|K\2=N > DS (-1)?*h}|K\{p+q)!/p!q! (B9) 
p ) q 0 
is a dominant series for (B8), and (B8) is uniformly convergent. It has 
been tacitly assumed that h + K 1. For other values of h and K, 
f must be analytically continued. ‘This question will not be discussed here. 
If g(t) = t, f(t) = t'(1+hm"+ Kn’)+. Thus, if g(t) may be approximated 
p 


by > a,t’, f(t) may be approximated by 
r= l 


S b,e[1+hm'+ Kap. 


i l 


If g(t) is an analytic function in /, f(t) is an analytic function in J. 


APPENDIX II]. ‘TABLE OF PERTURBATIONS 

rom the equations which govern the flow, many relations between the 

perturbed quantities may be obtained. A few of these are tabulated below, 

where all perturbations will be denoted by a bar, e.g. ¢, where the basic 
quantity is ¢. 

cy = yexp[(s, —s*)/c,]pi—, Pr = 22,4 /al(y—-1l—-piSe(y—-1). WD 








cy = yP\/p4; yP, = chp, + 2c, p, &. (11) 

’; = P2- i 2 Py _ P2— 91 Pi “Ps (111) 
Pi pi —9pz Py, Py py—9p, pi, — 9p» 

P,/P, =(14+ 9)|u, — W)2/e2? — 0. ({11’) 


P.+0P, P, ke -W ey 


awe, f. "laces 





P, P, aF. i OP,) uy 1. W Cy , 

=> -> = = -~—|]. \ 
“; -% P, E th 2] ti 
po mW a ss 
Pre" Pz t,—-W py uy -W" 


2(1 — 6)(W —u,)(W — ai,) — (u. — u,)(w — u,) x 
«| ae : cd =| = 2(1—0)c,é,. (VI) 





P = exp[(s—s*)/c,]p’, c2 = yP/p = yP°-’ exp|[(s—s*)/c,]. 
2é,/e, = (1—y2)P,/P, + 5,/yep- (VII) 
cz = Ovi +(1—4)c? = 605 +(1 Oye. 
(u, — W)(tig — W) + (u,— W)(a, — W) = 20(u, — W)(a, — W) + (1 — 8) 2c, é, 
- 20(u,— W)(a,— W)+(1—9)c, 2. (VIED) 
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In some problems, it is convenient to express all perturbed quantities 
in terms of one parameter, c.g. as in the problem of §4 where W was used 
as the parameter. Another convenient parameter, which leads to a concise 
representation of the perturbed quantities, is M = v\/c. From relation IIT’ 
above, 


Pp, P, 4yM, M, 


























P, P, ~ 2yMR—(—1)’ uF) 
Now Po P,/P, +0 a 2yM?2/(y + 1) 
p, 14+0P,/P,  14+6{2yM,/(y+1)-@ 
implies that 
Po PA 4M, % 0; 7 
me tae oe ~ eaeiiaalenweilee se X 
pP2 py M[2 + (y— 1) MZ] = (*) 
in view of the relation pyv. = p,v;. From c? = yP/p, 
2¢2/¢,— 26,/¢, = P,/P,—P,/P,— [P2/p2— p1/py], 
or 
Ge _ 1 _ Hy, 4 M, (1 + yMi) Xx 
a @  YR_G@Mi-y+1e+e- 1)M]2" om) 
From (III), 
2yM,/(y + 1)—6 = [—-0+2yM2/(y +1)f°, 
or M, My — 1—2yM2) = M, M,(2yM? — y +1). (X11) 
From the equation of state, 
exp|(s2- 1) C;| ie (Ps/Pi)(P1/p2)’, 
or -t Fe oy »| 2 7 P] 
Cy P, 1 Lee pd 
it follows that. __ 4yc,(y —1)(1-— M2)M, . 
5 Wy My TID yD —_ 
From M, = 75/¢., 
M, = (/C€z)[(0o/v2) — (€s/¢9)], (XIV) 
6, = ti, — W, o, = 4,-W. (XV) 


For a particular shock problem, the basic flow (i.e. all terms without 
bars) and the perturbed flow (subscript 1) are assumed known. If the 
problem is to be solved in terms of M,, equations (IX) to (XIII) yield 
P,, ps, &, My, 8. From (XIV), @, is obtained, and a, follows from (XV). 
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SUMMARY 


A small-scale experimental model has been constructed for the 
study of steady-state slow free convection of water in saturated 
sand. ‘The convection field is confined between coaxial cylinders 
maintained at constant temperature difference, and the upper 
and lower boundaries consist of planes with thermal and fluid 
insulation. Measurements of the temperature distribution 
within the convection space have been obtained for average 
boundary temperature differences (7';—7y) of (I) 18-40 A, 
(11) 32-70°A, and (III) 46°68 A. ‘Theoretical temperature 
values predicted from perturbation theory have been fitted by 
least squares. ‘lhe first-order estimates of 7/(7'— Ty), where 7 
is the convection parameter or modified Rayleigh number, do not 
differ significantly from a constant value (+3°, S.D.) for the 
three given values of 7, — 7), indicating good agreement between 
theoretical and experimental results. First-order estimates are 
made also of the temperature coefficient of thermal conductivity 5 
of the sand—water mixture, and of the coefficient of radiation loss ¢ 
at the upper insulated boundary, but these estimates are less 
reliable. Separate determinations of »/(7— 7 ), 6, c by direct 
physical measurement are in good agreement with the least-squares 


estimates. 


1. INTRODUCTION 
In a recent paper (Wooding 1957), the equations governing slow free 
convection of liquid in saturated homogeneous permeable material have 
been solved approximately, using perturbation expansions of the dependent 
variables (temperature and stream function) in powers of a convection 
parameter 7. ‘The method is very closely related to the technique of 

expansion in powers of the Rayleigh number 

x(T — Ty)gd® 
R = —————_ (1) 

Ko Vo 


* Now at the Cavendish Laboratory, Cambridge. 
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which has been applied to problems of free convection in fluids (Batchelor 
1954), and in fact 

ka( 7’, — T,)gd N(A/d)? 

ee | i) Prema) ay 


Under non-isothermal conditions each of these parameters arises naturally 





R. (2) 


when the corresponding equation of motion is reduced to dimensionless 
form—the Navier-Stokes equation in the case of R and the Darcy equation 
in the case of ». Hence 7 may be termed a modified Rayleigh 
number. 

In equation (1), « is the linear coefficient of thermal volume expansion 
of the liquid, 7, and 7’, are two reference temperatures on the absolute 
scale, g is the acceleration due to gravity, d is a representative linear 
dimension of the convection field, «) is the thermal diffusivity of the liquid 
at temperature 7), and v, is the kinematic viscosity at temperature 7). 
In equation (2), (K,,,/K,,.)y is the value at 7, of the ratio of the thermal 
conductivity of the permeable solid—liquid mixture to the thermal con- 
ductivity of the liquid; A is a linear dimension proportional to the particle 
or pore size of the solid, and N is a dimensionless constant dependent upon 
the geometrical shape of the particles. ‘lhe permeability k (= NA?) is 
assumed to be constant. 

Under given boundary conditions, the value of the convection parameter 7 
may be estimated in several ways. For instance: (a) by physical measure- 
ment of the quantities involved in the definition of 7 (equations (1) and (2)), 
(b) by measurement of the distortion of the temperature field caused by 
convective flow of liquid, the observations being compared with a 
perturbation model of the form 


fe 6, + 70, t n* 6, ates (3) 


0 
where @ is a temperature parameter and 4), 6,, 4, ... are coefficients calculated 
from perturbation theory, (¢) by measurement of the liquid flow rates at 
selected points of the convection field, the observations being compared 
with a corresponding perturbation model involving the stream function w. 
Methods (6) and (c) are applicable only over the range of small » for which 
perturbation theory is valid. 

This paper describes the measurement of 7» by methods (4) and (a) 
respectively, the former method using an experimental model with water- 
saturated sand as the convection medium, while for the latter method, 
separate measurements of the fluid permeability and thermal conductivity 
of the saturated sand are obtained. 

A particularly simple arrangement with cylindrical symmetry is chosen 
for the boundary conditions of the model. No attempt is made to use 
boundary conditions resembling those of a geothermal area, and, in fact, 
the latter task has been carried out recently by Mr John Elder in the 
Cavendish Laboratory, Cambridge. ‘lhe present model, in addition to 
providing a preliminary check of the perturbation theory, serves to indicate 
a practical method of measuring directly the ratio (K,,,/K,.)o%o Vo/R (see 








584 R. A. Wooding 


equation (2)) for samples of material obtained in geothermal areas, although 
the effect of large-scale fissuring, which is normally present in active parts 
of such areas, would be neglected by the method. 

‘Vo facilitate an accurate fit of the perturbation convection model to the 
experimental data, a refinement to the equation of energy transport, 
involving modification of the coefficient of thermal conductivity of the 
mixture, will be considered. ‘This is the variation with temperature of the 
thermal conductivity of the solid—liquid mixture. ‘The law of temperature 
dependence can be taken to be 

K,, = (Ky )ot 1 + (1— T 5, (4) 
where K,,,, (K,,,)y are the values of the thermal conductivity of the saturated 
permeable material at the temperatures 7 and 7) respectively, and the 
constant 6 is the linear coefficient of thermal conductivity. As 6 can be 
assumed to be small, its effect for small temperature differences can be 
treated by first-order perturbation theory. 

The diffusion or dispersion of heat due to percolation of the liquid about 
the sand grains can be expected to be a small effect, although it could become 
significant in cases of appreciably high fluid flow speeds about coarse solid 
particles. ‘lhe phenomenon closely resembles Brownian motion in that 
a given fluid particle performs small random motions which are superimposed 
upon the main fluid motion q. From the analogy with the classical theory 
of random flights (Rayleigh 1899, Chandrasekhar 1943), it is a simple matter 
to show that the magnitude D of the ‘thermal dispersion coefficient’ will 
obey the relation 

D = O(pcd q)), (5) 
where A is a typical dimension of the solid particles or of the pores, pc is 
the heat capacity of the fluid, and q, is the main flow speed. ‘Then the 
apparent thermal conductivity of the saturated permeable material, for 
isotropic thermal dispersion, becomes A,,+D. If the dispersion eftect 
is non-isotropic, its coefhcient assumes a complicated tensor form 
(Chandrasekhar 1943). However, for very small A in comparison with 
the boundary dimensions, and with very slow fluid flow, as in the present 
experiment, it can be shown from (5) that D is negligible in comparison 
with A,, and the etfect of dispersion can therefore be neglected. 

For purposes of the experiment, it is necessary to consider briefly the 
influence of the dynamic boundary layer associated with each rigid boundary 
which is immersed in permeable solid material. Darcy-type flow solutions 
cannot satisfy the ‘ non-slip’ condition at a rigid boundary, whence it becomes 
necessary to neglect the region of boundary layer flow. ‘his approximation 
is justified when the Reynolds number of the flow is very small. ‘Uhen 
all fluid velocities in the pore spaces are of laminar (Poiseuille) type, and 
the shape of each velocity profile across a pore space, or normal to a rigid 
boundary, is independent of the kinematic viscosity v and of the main fluid 
How rate q. It follows that the typical thickness 6 of the boundary layer 
depends on only the typical solid particle (or pore) size A, and that the ratio 
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6/\ is of order unity. ‘This follows from the particular form of the 
Navier-Stokes equations for which the inertia terms are negligible, gradients 
perpendicular to the boundary are large compared with those parallel, and 
the pressure gradient is given by Darcy’s law. If X is a coordinate measured 
along an impermeable rigid boundary immersed in saturated permeable 
material and Y is measured along the normal, one has 


, 


(pressure gradient)/p = — Z 4 = “yy 


where u is the X-component of the flow, and q, is the speed of the main 
stream outside the boundary layer. Since uw = O(\q)/6"), this gives 


5 = O(k"2) = O(A), (6) 
from the relation k = NA? (Slichter 1899). 


2. APPROXIMATE SOLUTION OF DIFFERENTIAL EQUATIONS 
As in a previous paper (Wooding 1957), the differential equations of 
slow steady flow may be expressed in dimensionless units. ‘The temperature 
parameter @ is defined by 
6 = (T—T))/(1,— To), (7) 
where 7,, 7, are two typical boundary temperatures. ‘The ‘equivalent 
flow’ vector q, is defined froin the continuity relation 
PF = Po Go» (8) 
in which it will be assumed that the density p of the liquid is related to the 
temperature by the approximate expression 
mi a ee ; 
P= potl — «(7 — Ty)? — B(T, — Ty) OS, (9) 


py being the density at temperature 7, and «, 6 being constants. ‘lhe 
modified Péclet number € is given by 

(KL /KujoS = Gp d/o; (10) 
where «, is the thermal diffusivity of the liquid at temperature 7), d is a 
typical length derived from the boundary dimensions, and (K,,,/K,,)o 1s, 
as before, the ratio at 7, of the thermal conductivity of the mixture to the 
thermal conductivity of the liquid. 

Since the boundary conditions of the experiment have cylindrical 
symmetry, it is convenient to take axes in dimensionless units 7, 3, with 
Or horizontal and Os directed vertically upwards, and to define Stokes’ 
stream function % by 

C = id,/r—ku,/r, (11) 
i, k being vectors in the directions Or, Oz respectively. ‘The suffices in r, 3 
signify partial differentiation. 

Then, for slow steady motion of water in a saturated homogeneous 
permeable solid, the equations of energy and motion give in cylindrical 


F.M. 2P 
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coordinates a pair of simultaneous partial differential equations for 6(7, 3) 
and d(r, 3), 


1 

6.4 1 fee eee b(T, = ry} (76..).. + (rG6.) = — (a5 d.—uw. 6.) 

/ r / om Y L / . - Y ~ / Tr ~ 
= Lipa (12) 
ro(r, 2) 

“2 1, } Gis 40) Ok +0. wb) 
aay l+a(7T,-—T,)@> ' 
nr{1+a(T,— 7,01 +(28/2)(7, — 7)0}8, = 9, (13) 
where 


a = (T, — 260)" ( A)! (14) 


is a temperature coefhcient for viscosity. 

For the approximate solution of equations (12) and (13) it will be assumed 
that », b are small, and that #6 and y may be expanded in power series 
involving these quantities. (Also, for convenience in later analysis, we 
introduce the notation m = y/(7,— 7), the parameter m being independent 
of T,—T,.) Then the perturbation models will be written as 


6 = 0,4+m(T,— T,)0,+6(T, — T))0 + mT, — Ty)*0. +... (15) 
j= att Tb Oh, OP... (16) 


Here, the second-order term (in m”) will be calculated in order to estimate 
the magnitude of the second-order correction. ‘The second-order terms 
in mb, 6? will be ignored, on the assumption that only the linear term in 6 
has significance. ‘The series solution (15) and (16) of the equations (12) and 
(13) may be fitted to suitable experimental data to obtain estimates of the 
parameters m and 6. 

Substitution of equations (15) and (16) into (12) and (13) gives a series of 
partial differential equations for the perturbation coefficients 6, %,, 44, 


bo, 0, and 6%, 


(Laer = ade + Padee— 110) Pa) + (00) a) e3 
t ril r a(T, T)%} ul u (2p o)( rT, ” T)9y§(90) + = U0, (18) 





1 A(8p, 4) 
AO) + 218) 2 oe 9) 
(73) wr 1) r ar, 2) (19) 
1 
(tbo)-.— — (Yo) (eb) YF) (Yo) + (Fo) (%).} — 


= YF) (4a) + (PrdCbrde§ + Y7F (Fo) (41), + (Go)e(r)-5 4 

tr{1 +a(T,— Ty)6ot{1 + (2B/a)(T, — Ty) 65 }(9,),.4 

+ rf{1 + a(T, — T)9o}(28/«)(T, — To), 4 
+{1+(2B/a)(T, — T)0}a( 7, — T)91\(8o)» = 9, (20) 
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where y = a(T,— T,){1 + a(T,— Ty) 45, 

l 1 O(8, a) — 10(,, By) 
Os (IO a See eee 21 
(Vo) r, yr’ 2) ( 2) 22 r a(r, z) r O(r, 2) ( ) 
y(1) I Aa )(1) §(A \2 A)? } JA) 
Or + Or +e or + (odes = 9. (42) 
/ 

‘These are readily transformed into finite-difference equations by the usual 
procedures, taking a square mesh of element dimension Af, and relaxation 
techniques may be applied in order to obtain approximate numerical 


solutions for the coefficients. 


3. DESCRIPTION OF APPARATUS AND EXPERIMENTAL PROCEDURE 
Figure 1 illustrates in diagrammatic form the arrangement of temperature 
and flow boundaries chosen for the experimental work. Oz is the axis of 
two concentric cylinders 4D, BC which have radii in the ratio 2:7. In 


y 


A 


i] 





z=08 











the cylindrical model, showing mesh points of 


Figure 1. Boundary configuration of 
the relaxation net 


the dimensionless system the length unit is taken as d= AB, giving the 
dimensionless 7, z values shown for the various boundaries. ‘The 30 points 
arranged in a square mesh represent points of the relaxation net, from which 
it will be seen that A = } in dimensionless units. ‘These points correspond 


to thermocouple positions in the experimental arrangement. 
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The cylinders AD and BC are held at constant temperatures 7, and 7) 
respectively (7; T,), while the boundaries AB, CD satisfy conditions 
of nearly complete thermal insulation. All four boundaries are impermeable 
to liquid. Consequently, the boundary conditions for the numerical 


calculations may be described ideally as follows: 


on BC, G= 9; 

on DA, d6=1; 

on AB, CD, ¢, =O; (23) 

on AB, BC, CD, DA, = 0. (24) 
\ feature of these boundary conditions is that (4), vanishes throughout 
equations (17) to (22). 
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i 2. Leneth-breadth distribution of 100 grains of quartz sand. ‘The regression 
lines and standard deviation have been calculated from a total sample of 233 
1 the experimental arrangement, two copper cylinders of radii 3-8 cm 
13-3 em are arranged coaxially to represent boundaries DA and BC 
spectively. ‘lhe outer cylinder is surrounded by a water jacket, while 
the inner cylinder forms re pr of a water-tight container. Constant, 


a9 . : 1 aa a ake 
.ctically uniform, temperatures are maintained on the boundary cylinders 
yy means of circulating water with thermostatic control. 


The space between the temperature-controlled copper cylinders, of 


9-5 cm, is filled with standard graded quartz sand (porosity 0-4 
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approximately) to a depth of 0-8d. Figure 2 is a scatter diagram showing 
the length—breadth relation for a sample of 100 sand grains, based on the 
assumption that each grain is approximately ellipsoidal. (Since these results 
have been obtained from a single photomicrograph, measurements are not 
available for the depth of each grain.) ‘The sample is seen to be satisfactorily 
homogeneous, and if A is a typical dimension for grain size it follows that 
A/d = O(10-?). 

A thin sheet of bakelite mounted vertically and radially supports 
copper—constantan thermocouples in the 30 positions shown in figure 1, 
the thermocouple leads being carried circumferentially for several centimetres 
from the junctions to reduce conduction errors. Each thermocoupl 
junction is arranged to project about 0-4+cm beyond the bakelite sheet, 
and is supported by gluing the copper and constantan wires to the sheet 
For the thermocouples adjacent to the boundaries corresponding to B 
and )A, extra support can be afforded by the cylinder walls, as these are 
in effective thermal contact with the thermocouples. 

\lthough a boundary layer exists close to the vertical bakelite sheet 


y } ] . ] ’ } +4 - } sank 
its effect upon the thermocouple readings can be considered negligible 
l 


provided that 6 == A 0-4 cm, where 6 1s the thickness of the layer (§ 1). 


From figure 2 it appears that this condition is practically satisfied. 


Thin sheets of bakelite sealed with a rubber-based glue render the upper 
and lower boundaries impermeable to moisture, and are backed by layers 
3 cm thick oi P< rspex ‘honeycomb’ to pros ide thermal insulation. 
Compression pphi tO ton 1@ upper lavet ( nsulating material 
is tran sughout t] sembly, and ensures t except for possibl 
packi: regula pA inc oul in approximatel ist 
Huid pei bilit | thermal conduc ty within the convectio yACE 

\fter the convection space has been filled with boiled water, care being 
taken to ex lude all air bubbles, the inner and outer boundaries are adjusted 
to h: suitable temperature LitT¢ 1 which j intained s nearly 
constant as possible \ pr id of up to four hour elapses perore steady-state 
conditions are approached 

hree different values of the boundary temperature difference have beet 
employed, the experimé tal average values of 7, T, be ing (1) 18-40° A, 
(II) 32-70 A, (III) 46-68 A. In each case, 7', is held close to (273 + 20)° A, 
for which « = 2:25 x10? (° Ay’, 6 3-75 x 10-* (© A) in equation (9), 


and a = 1/33(° A)! in equation (14). 

Experimental values @,. of the convection temperature parameter are 
computed from the measured temperatures 7'( A) using the formula (7). 
The results for the three experiments are given in table 1. 

Steady-state conditions are found to be attainable on each cylindrical 
boundary to within +0-05 \ for experiments | and Il, and to within 

():20 A for experiment III. ‘These results are obtained by observation 
of several boundary thermocouple outputs whilst the thermostats controlling 
the boundary temperatures pass through several ‘on-off’ cycles. It ts 
easily shown from formula (7) that resultant variations in time of the 
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parameter @, are less than + 0-005 for experiments I and II, and less than 
+ 0-010 for experiment II]. Corresponding time variations in the interior 
thermocouples are much reduced, owing to the short period (five to eight 
minutes) of the boundary variations, which are lagged by the saturated sand. 








2D 999 C686. AG PAZ? «Ay O7s: <2) te 3) ae tf “th 
Expt. 986 =(14) | 681 (4) | 459 (2) | 279 (1) } 119 (1) QO (0) | 
I | £ 1010 (—10) | 673 (—1) | 445 (0) | 263) «~(0) | 112 (0) | F 1 (—1) 
| 1004 (—4) | 655(—1) | 416 (0) | 243 (0) | 100 (—1) 1 (—1) 
1 998 (2) | 645 (0) | 399 (O) | 241 (0) | 103 (—2) | B 6 (—6) 
| | ae eS eee . ae 
| D1003 (—3) | 735 (0) | 514 (0) | sis “<(h) | 148 (2); C 6 (6) 
Exot. | 997 (3) | 710 (0) | 486 (0) | 305 (0) | 141 (1) ; (2) 
7 k 1009 (—9) | 680 (—2) | 440 (—1) | 263 (—1) 124 (—2) | F 7 (—7) 
999 Ch) | 613 <0) | 387 (©) | 224 (—1) | TOL (1) t 
A 993 (7) | 587 (2) | 353 (0): | 210\(—1) | 98 (—1) | B 2 (—2) 
| D1020 (—20) | 791 (—8) | 593 (4) | 368 (3) | 189 (—3) | ¢ 3 (—3) 
: | 1008 (—8) | 742 (—5) | 533 (—3) | 293 (—3) | 166 (—3) 5 (—5) 
Expt. | ,, - ri Sa aes ‘ ' é 
II £1009 (—9) 600 ( 5) | 437 (—1) | 262 ( ) | 123 (—2) | J o | 3) 
; 99] (9) | 559 (4) | 347) (2) | 197 (2)} 86 = (3) G2) 
1 973 (27) | 509 (10) | 295 (4) | 172 (3)| 76 (5)| B-—10 (10) 
Table 1. Experimental values of the temperature parameter 10°@,. Corresponding 


values of the boundary correction term 10°@) are shown in parentheses. 


However, in the steady state, small spatial temperature variations occur 
in the Oz-direction along the boundaries DA and BC (table 1). ‘These 
appear to be due to the finite effective conductivity of each cylindrical 
boundary system, emphasized perhaps by the presence of a steady flow 
pattern within the circulating water. 


4, ANALYSIS OF RESULTS 
As noted at the end of §3, the steady-state values of @, obtaining on 
boundaries BC and DA vary by small amounts from the theoretical values 
of @ given by (23). The theoretical boundary conditions will now be fitted 
to the experimental values by applying a perturbation 6) to 6, such that 
on BC, DA, 1—@© = ae 
on AB, CD, 6 = 0. (25) 
Also, let the definition of 6 be such that 
6-0 = 6, (26) 


within the convection space. Here, @ can be assumed to be represented by 
the perturbation series (15). As 0 is very small, cross-product terms 
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involving m(7', — Ty), 0, etc., will be neglected, whence it is readily shown 
that 6 is harmonic to first order in the convection space. Values of 6 
calculated by relaxation methods are given in parentheses in table 1. 

A further correction to the boundary conditions arises from the possibility 
of thermal radiation through the insulation at the boundaries 4B, CD. 
Since the external ambient temperature is normally close to 79, it is possible 
to assume that the radiation boundary conditions take the form 

on. CD, 0. = co: 
on AB, 6, = c’é. (27) 
With the assumption that the radiation constants c, c’ are very small, it is 
possible to take first-order perturbations in @ of form 
cB + ¢'G®), (28) 


this expression being combined with (15). 








| 10°, | D 1000 | 677 447 | 269 123 ¢ 0 
a ee eee | =) : 
10%, | D 0 | oO | 0 | Oo | 0 ( 0 
| 0 | 41. | 52 48 | 32 0 
| | E 0 | 51 | 67 | 62 | 42 | 1 0 
| | | | 
1049, | D 0 131 | 136 | 96 7 ( 0 | 
0 90 94 67 ao Tl 0 
| fi: 0 0 0 0 Ooty 0 
10°ubso | D 0 | 0 0 0 0) ( 0 
| O 8 30 46 38 | 0 | 
/ 0 0 | 0 | 0 0 |Z 0 
10°, D 0 25 2 17 | 13 ( 0 
| 0 | 27. | 5 1 | 10 0 
| J 0 28 | i || 5 7 F 0 | 
} | 
ra, oe parr wane = [ | 
| 10°90) D 0 | 112 | 126 | 100 55 ( 0 
ee | ee eens = | “ = = | [- = 
| 10390) =| p 0 | 144 | 146 | 105 | 52 | ¢ 0 
0 | 68 sO. 62 32 0 
} 7 | E 0 | 35 45 | a 20 | F 0 | 
| | 0 | 21 28 24 | 13 0 
| | 0 | 17 23 20 | 11 B 0 


| 
| 
| 
| | 
| | 





Table 2. Numerical values of the perturbation coefficients calculated by relaxation 
methods for experiment I. The given values of 6), 6“) and 6) apply for 
experiments I, II and ITI. 


Tables 2, 3 and 4 give numerical values of the coefficients 4p, #,, 0, who, Ao 
and 6, calculated by relaxation methods from equations (17) to (22) 
respectively, and using the boundary conditions derived by substitution 
of (15) into (23) and (16) into (24). Also tabulated are values of 0°, which 
are calculated from Laplace’s equation with boundary conditions derived 
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by substituting (28) into (27). In all cases the relaxation net of figure 1 
is used. In general, it is necessary to tabulate only part of the calculated 
field for each case, as both @, and 6 are constant with respect to z, while 
%, and 6, are symmetrical about the horizontal plane EF’, and 6, and x, 





| | 
10%, D 0 0 0 0 | 0 | Cc 0 } 
| 0 59 CI 70 | 60 | Ss 0 | 
a | na a ie 
I 0 74 | 91 79 | 50 | 0 | 
= : | ; . : eee | —|—___—— 

1040, D 0 184 184 127 | a ic 0 | 
0 | 126 127. | oe | 43 | 0 | 

I 0 0 0 0 0 | 0 

| | 

10° u/s. D 0 |} 0 0 | 0 | 0 ee 0 

0 0 77 93 65 0 

I 0) 0 0 0 0 I O 

10°54 D i) 44 9 33 | 24 | ( 0 

0 49 4 | 2 19 0 

E 0) 53 16 ign i 14 |F 0 


] 


Table. 3. Numerical values of the perturbation coefficients calculated by relaxation 


P OS acraiae ae | 
methods for experiment 11 


LO () Q) 0 Q) 0) ( 0 
0 S( () 74 $4 Q 

/ () (0) 117 Q7 38 } () 

()+f ie Q) +> 2 38 160 / ( () 
) L6 164 111 5 0 

} () 0 0 0 Q i {) 

10 D () () () () 0) (” 0 
) 2? 158 165 101 0 

i 0 0 0 0 0 F Q 

L0°¢ D 0 70) 23 58 39 (" 0 
() Pel uv 4( 5 QO 

} O 90) 22 23 23 I @) 


Table 4. Numerical values of the perturbation coefficients calculated by relaxatiot 
methods for experiment ITT. 


are antisymmetrical about EF. 6° does not possess symmetry of this type, 
but 6° and 6®@ are related antisymmetrically about EF. Finally, the 
functions 4, 6 and 6 have the same values for experiments I, II and III. 

Because of the small scale of the apparatus used, the sand grains are of 
appreciable size (A d = O(10°*)) and the spatial temperature gradient is 
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high, with resultant scatter in the experimental temperature readings. It is 
preferable, therefore, to use statistical methods (least squares) in fitting the 
theoretical model to the data of table 1 for the three values of the boundary 
temperature difference which are employed, and to use significance tests 
in evaluating the adequacy of the model. 

Combination of equations (15) and (28) shows that the statistical 


temperature model will be of the form 
er, mT, — 7,)6,+4(7, T,,)O + cB + c'@® 4 mT, Ty)?4.+€, (29) 


to the order of accuracy considered, where all of the parameters m, 4, « 


and ¢’ should be constant for varying values of 7— 7). Point values of 


the residual error variable € will be assumed to be distributed approximately 
normally and independently of each other with zero mean. Since not all 
terms in the analysis are necessarily significant, and since the complete 
model would prove unwieldy when obtaining a least-squares fit, a series 
of partial models will be used, each model involving not more than three 
of the parameters on the right-hand side of (29). For future reference, 
each partial model will be labelled H;, with residual variable e; and 
. ; 


corresponding sufthix will be used to identify the parameters retained. 


1° = ra | P aa a +14 . smart eh « shar . 

lable 5 gives the numerical results obtained with a number of partial 

j 1 . 1 1 = +o i 1 ] rt kk Lia 4 . . 1 ] 
models for the three cases I, 11 and II] bv the 1 | least-squares methods 
In table >, entries from the left are in terms of model numbe1 7, Darameters 
retained 1n model 7, degrees of freedom v,, ex eriment number ({, LI o1 iil 
optimal parametric estimates, and standard error of estimates 
: oe cg eonil id Diane repay catia ¢ } . . i : 
signincance of eacn added paramere! is estimated DY iCallS ¢ the variance 

. ’ 1 . y 1 “1 
ratio /’... computed from the residual sums of squares >(e7), +(e7) for the zt] 
. i a tot 

ne ) 1 ? 1 } 

and jth models respectively (2 > /), from the formulk 
iN (e2) (e4)'/( ) 
fk KO y:) = — : {rf 
} , 
N(e2) 

r) 1 400 A4\ } 1 | sera ed 
(Pearson & Hartley 1954). Each tabulated value of O; gives the probability 


of exceeding the observed F-ratio by chance. ‘Thus, a very small value 
of O implies that the added parameter 1s significant. 

lhe models H,, H, and Hy, (2 1,2, 3 in table 5), are intended to test the 
significance of the convection effect to first and second orders respectively. 
However, H., cannot be tested against H, by the usual significance tests, 
since m, m® are not independent parameters. A conservative test for 
non-significance of m* may be obtained by replacing m? with an independent 
parameter /(H,), and testing H, against H,. ‘Thus, if this test shows that 
the effect of / is non-significant, then the effect of m? will be non-significant 
also. 

Values of the parametric estimates *, = jj/(7,— Ty), @5, bg from table 5 
are plotted in figure 3 (a), (6) and (c) respectively, these being chosen since 
they represent values estimated from the simplest models. 
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Also shown in figure 3 are overall estimates mg, ¢;, ., which are obtained 
by assuming the hypothesis that the quantities m, c, 6 are constant for 
experiments I, If and III, and obtaining least-squares estimates using the 
‘These results are summarized in table 6, and are dis- 


combined data. 
The variance ratio F° 


= tinguished from the results of table 5 by dashes. 
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Figure 3. Values of the parametric estimates from table 5 for experiments I, I] 
and III, and the overall estimates froi.. table 6 for (a) m7, () ¢5, (¢) 6. 


is used to examine the significance of the effect of each added parameter 
m., b., c., and the variance ratio F;; is used to examine the significance of 
the increase in the residual sum of squares }(e;2) over the sum of the three 
individual sums of squares ¥(e?) given for the 7th model in table 5. 
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5. DISCUSSION OF STATISTICAL RESULTS 

In table 5 for ¢= 1,2, the values of F, show that the parameter 
m = y/(T,— 7.) in the first-order model has a very high significance, which 
increases with (7',;—7)) as the convection effect becomes increasingly 
dominant over secondary effects. Hence the null hypothesis H,(4— 4, = €,) 
must be rejected, and the presence of a convection phenomenon appears 
to be well established. However, from Fy, the effect of the second-order 
term (m?) is not significant in experiments I and II, and the effect in III 
is probably not significant. Presumably, the term in m? would become 
significant for larger values of (7',—7,) than are available in these 
experimental data. 

Separate tests of the effects of the parameters m, c and m, 6 against m 
(F'52, Fy. respectively) show comparable significance for the reduction in 
residual sums of squares in the two cases. However, with the ettect of m, 
c, 6 tested against m, b (F;,), there is not a significant reduction in residual 
sum of squares. It appears that this result is due to the similarity, analogous 
to a high correlation, between the coefficients of c and 6. In physical terms, 
the effect of radiation through the upper insulated boundary (coefficient c) 
is very similar in form to the effect of variation in thermal conductivity 
due to the negative temperature coefficient (b) of the quartz sand. Both 
phenomena tend to reduce the temperature within the convection field, 
and since each effect is small, the present method of analysis is not able 
to resolve them completely. 

An additional test, of the parameters m, c, c’ against m, c, shows that c’ 
has no significance. It follows that radiation through the lower boundary 
surface is negligible. 

‘To examine the significance of the dispersion term (§1), a further test 
is formulated using m, 6, ¢ and an additional perturbation parameter 
representing the effect of simple isotropic dispersion. As the three values 
obtained for the coefficient of the additional parameter are scattered closely 
about zero, it can be concluded that there is negligible dispersion of heat 
by irregularities of the fluid motion. 

Considerable heterogeneity of residual variance d(é v;) has been tound 
between the results of experiments 1, II and III, much of the heterogeneity 
being due to a relatively large residual variance in II]. Physically, it would 
appear that the higher temperatures involved in III have resulted in poorer 
stabilization, the residual standard error in the temperature measurements 
being given by (1) +0-12 A, (I1) +0-18°A, (III) +0-69° A respectively, 
after m, 6 and c have been fitted. ‘hese values are larger than the time 
variations mentioned in §3, since steady-state scatter in the thermocouple 
readings is present also. ‘The effect of this heterogeneity of variance is to 
modify the distribution of the F’-values given in table 6 and the corre- 
sponding values of O,, may therefore be only an indication of the true values. 
Although it follows that the results of table 6 are statistically less reliable 
than the results of table 5, the Q,-values in table 6 are all either so large or 
so small that any numerical modification necessary is unlikely to alter their 
physical implications. 
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The results of the overall tests F3,, F'5;, Fy, 5, (table 6) show that the 
hypothesis that m, 6, c are constant is not disproved. ‘lhe standard error 
of the estimate m, is about + 3°... Constant values of these three parameters 
are expected, of course, since none is dependent upon (7', — 74), the quantity 
varied between the three experiments. 

It will be seen that the test F',,, (which tests m, b, c against m, 6) does not 
give significance. If the simpler model is desired, therefore, one could 
ignore boundary radiation and choose the model H,, in preference to H>. 

A further hypothesis H,, using the overall data to fit m to second order 
and 6 to first order as in equation (15), is found to give the same residual 
sum of squares as H,. Hence there is no preference for or against the 
second-order model when it is compared with the first-order model. 


6. SUPPLEMENTARY OBSERVATIONS AND RESULTS 
In order to calculate an approximate value for m by alternative means, 
it is necessary to measure the permeability of quartz sand under compression 
conditions resembling those of the convection experiment, and to measure 
the thermal conductivity of quartz sand at temperature 7 in the absence 
of convection. A standard permeability experiment of Darcy type (Muskat 
1937) has been set up, the sand being carefully tapped down and compressed 
as in the convection model, and a value of 

k = (261 +005) x 10° cm? (31) 
has been obtained. For the conductivity measurement, water-saturated 
quartz sand is compressed between a pair of square horizontal plates, the 

upper plate being at the higher temperature. A conductivity value of 
(K,,, )o = (4:60 + 0-05) x 10>? c.g.s. units (32) 


Jo» and 
1 


has been obtained at (20 +273) A. With these values for k and (K 
with « = 2°25 x 10-*(°A)“!,.d = 9-5 cm, g = 980. cmsec™*, ») = 0-010 cm* sec” 
units, equation (2) gives 


m = n/(T,—T) = 0-119 + 0-003 (° A). (33) 


m 


(‘his value is in very good agreement with the least-squares estimates of m 
in $4, table 6. It will be seen that m has been determined with approximately 
the same standard error by either method. 

l'rom these values of m, it is a simple calculation to show that the ratio 


(K,, K,,,.)o(Ko ¥o/R) = 17:5 + 0°5 cm? sec (34) 


for 7, = (2734+20)° A and tor the particular water-saturated quartz sand 
used in the experimental work. ‘This result, which is a property of the sand 
and water alone, is of more general usefulness than m, as it is independent 
of the boundary conditions used in the model. 

A rough estimate of 6 has been obtained from the data on thermal 


conductivity of quartz published in the International Critical ‘Tables 
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(National Research Council, U.S.A. 1929). In the temperature range 
20°C to 70°C, the thermal conductivity of quartz can be represented 
approximately by 

K, = 0-029/1 — 0-0037(7T— Ty)! c.g.s. units 
parallel to the principal axis of the crystal, and 

K, = 0-016/1 —0-0024(7'— T,)} c.g.s. units 
perpendicular to the principal axis, where 7) = (273+20)°A. If the 
directions of the principal axes of the sand grains can be assumed to have 
. spherical distribution, and if the quartz contributes the major share of 
the thermal conductivity, the average value of the temperature coefficient 
would be approximately 

b 0-0028(° A)-}. (35) 

This is to the same order as the values of 6 estimated in § 4, table 6. 


1” 
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‘Table 7. Measurement of flux of heat output through upper insulation. 


\n independent estimate of c can be obtained from the approximate 


formula 


p a ~]1-4 “1-4 
2rd 
ns er ae ae (36) 
(K,,Jo(11- ly) O-4 J O-4 
where g is the outflux of heat per unit area in c.g.s. units, and 6® is the 


coefficient of c in the perturbation formulae. ‘lo measure q, a series of four 
differential thermocouples (calibrated as heat flux meters) are placed on 
the upper insulation at radii r corresponding to points of the relaxation net. 
The results are given in table 7, using the boundary conditions as for 
experiment I. After the integrals in (36) have been evaluated approximately, 
a value of 

0-10 + 0-05 (37) 


is obtained. ‘This result is to the same order as the values of c estimated 
in $4, table 6. 
7. CONCLUSIONS 
Evidently the first-order perturbation scheme provides a satisfactory 
fit to the convection data for the range of boundary temperature differences 
(7,—7 ) chosen. Its range of validity extends from 7» = m(7,—7)) = 0 
to the point where the second-order term becomes significant in comparison 
with the term of first order. Presumably terms of higher order can be 
calculated, with increasing computational difficulty, but accuracy is not 
easily maintained owing to cumulative errors. 
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With the above limitation to first-order accuracy, it appears that the 
following additional conclusions have been justified experimentally. 

(a) The motion of the liquid percolating through the sand under 
non-isothermal conditions is adequately described by Darcy’s law for the 
low values of Reynolds number (< 1) applying here. 

(5) ‘The hypothesis that the packed sand behaves as an approximately 
homogeneous permeable medium is satisfactory. 

(c) Near a fixed impermeable boundary, changes in sand packing 
due to the presence of a rigid surface might be expected to produce a thin 
layer of altered permeability. Neglect of this effect does not appear to 
introduce detectable errors: hence it is not of importance. 

(d) Modifications to the apparent thermal conductivity at very low 
Reynolds numbers by fluid dispersion (§ 1) is negligible. 

(e) The thickness of the fluid boundary layer adjacent to smooth rigid 
boundaries is small, probably close to the typical grain size. 

‘The conditions (a) to (e) represent initial assumptions which are 
contained implicitly in the theoretical solutions. Departures from these 
conditions are negligible in comparison with the convection and other effects 
analysed above, since the residual variance after fitting m, 6, c is very small. 

While the analysis of thermal boundary loss and variation of thermal 
conductivity with temperature does not give accurate results, it is clear 
that these effects are small when compared with the convection effect. 
he estimate of m is not altered appreciably by their presence, and the 
general conclusions about the convection process are not affected. 


-~ 
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SUMMARY 


\n estimate is given of the distribution of skin frictional force 
per unit length, and of displacement area, on the outside of a 
semi-infinite cylinder, of arbitrary cross-section, moving steadily 
in a direction parallel to its generators. A Pohlhausen method 
is employed with a velocity distribution chosen to yield zero 
viscous retarding force on the boundary layer approximations. 
(The smallness of the fluid acceleration far from the leading edge 
has been pointed out by Batchelor (1954).) Like the Rayleigh 
method, this method is expected to yield reasonable results at 
large distances from the leading edge. However, for a large 
class of cross-sections, which includes all convex cross-sections 
and locally concave cross-sections with re-entrant angles greater 
than }, the method yields the expected square root growth of 
the boundary layer at the leading-edge, with a fairly close approxi- 
mation to the coefficient, and it is supposed that the skin-frictional 
force and displacement area are given with reasonable accuracy 
along the whole length of the cylinder. 

Results for the elliptic cylinder and the finite flat plate are 
given in closed form, valid for the whole length of the cylinder, 
and are expected to be in error by at most 20°... In addition, some 
estimate is given of the effect of corners on skin frictional force and 


displacement area. 


1. INTRODUCTION 

This paper deals with the boundary layer on a long cylinder of arbitrary 
cross-section in a stream parallel to the generators. ‘he shape of the nose 
is immaterial, but it is supposed to be smooth enough not to cause separation. 
It is the growth of the boundary with distance x from the leading edge, and 
the distribution of skin friction along the cylinder, with which this papet 
is concerned. 

‘The problem for the cylinder of arbitrary cross-section has been tackled 
by Batchelor (1954). His work is based on Rayleigh’s (1911) method of 
inferring a rough answer from the analogous time-dependent heat-conduc- 
tion problem of determining the temperature distribution in an infinite 
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homogeneous solid, initially at zero temperature, but with the internal 
boundary (an infinitely long cylinder) maintained for t > 0 at a constant 
non-zero temperature. ‘The time ¢ is then interpreted as x/U, the thermal 
diffusivity as v, and the temperature as the departure of the velocity u 
from the free stream value U. 

This method assumes that vorticity in the boundary layer is convected 
downstream with the free stream velocity instead of the local velocity u 
in the boundary layer. In consequence, the boundary layer thickness is 
underestimated and the skin friction is overestimated. 

For the cylinder of arbitrary cross-section Batchelor finds that, for 
small values of vx/U, the local properties in the boundary layer are 
independent of the cylinder cross-section. They are as in the well-known 
Blasius layer, for which the boundary layer thickness varies like (vx/U)}*. 
For the Blasius layer, Batchelor, like Rayleigh (1911), obtains for the local 
skin friction 0-5644U(vx/U)~!*, when the true coefficient is 0-332. This 
is a serious overestimate. 

At values of vx/U for which the thickness of the boundary layer is 
comparable with a typical dimension / (taken to be the perimeter divided 
by 27) of the cylinder cross-section, the effect of the cross-sectional shape 
is important. ‘The first approximation for the effect of the shape of the 
cylinder on the force per unit length of cylinder is determined in terms 
of the number of corners and their angles in the cylinder cross-section; 
if there are no corners, the force on a unit length of the cylinder is the same, 
to this approximation, as that on a circular cylinder of the same perimeter. 
For large values of x/U the boundary layer thickness is large compared 
with /, and it is here that we expect the Rayleigh method to furnish reasonable 
results since conditions in a large part of the outer boundary layer approxi- 
mate to free stream conditions. Batchelor shows that the bulk properties 
of the boundary layer for large values of x/U are the same as those for the 
circular cylinder which has the ‘equivalent radius’ c, such that the two 
cylinders, if given the same charge per unit length, would have the same 
electrostatic potential at large distances. Moreover, the rate of retardation 
of the fluid becomes very small. 

Batchelor’s solution for the general cylinder can be regarded as giving 
useful qualitative results, but is inadequate quantitatively except far down- 
stream where it was shown by Glauert & Lighthill (1954) to be good. 

Glauert & Lighthill (1954) investigated the boundary layer on a circular 
cylinder by using a Pohlhausen method with a velocity profile of the form 
A(x)log(r/c), where r is the distance from the axis of the circular cylinder. 
This profile has the virtue of satisfying conditions near the wall as accurately 
as possible, and of being the profile for which the viscous retarding force 
is identically zero. A(x) was determined by substitution in the integrated 


momentum equation. 

‘The results near the leading edge show good agreement with the exact 
solution (due to Seban & Bond (1951) and corrected by Kelly (1954)), 
in which the Stokes stream function is expanded in powers of {vx/(Ua?)}?, 
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‘The solution also agrees well with an asymptotic series solution, published 
in the same paper, which is valid far downstream. 

To summarize: Batchelor has shown that for a cylinder of arbitrary 
cross-section far from the leading edge the viscous retarding force is very 
small, and in the outer part of the boundary layer where conditions are 
symmetrical the velocity profile is of the form A(x)log(r/c), where A(x) 
depends upon the cylinder cross-section. Moreover Glauert & Lighthill 
have shown that for the case of a circular cylinder a Pohlhausen method 
with a velocity profile yielding zero viscous retardation gives good results 
along the whole length of the cylinder. 

In the present paper a Pohlhausen method is again used with a velocity 
profile for which the viscous retarding force is identically zero along the 
whole length of the cylinder, and which agrees with the known solution 
far downstream. We expect the present method to furnish good results 
in the region where the boundary layer thickness is comparable with /, 
since in this region Batchelor shows that the viscous retarding force is small. 
‘To see whether the method, like that of Glauert & Lighthill, has value also 
near the nose, we expand the solutions for different shaped cross-sections 
at the nose and compare them with the expected Blasius solutions. 

We assume the velocity profile is of the form A(x)d(y, 3), where y, 2 are 
Cartesian coordinates in the cross-sectional plane. It is further assumed 
that for a large distance r from an axis fixed in the cylinder boundary ¢ varies 
like log(r/c). The determination of ¢ is then seen to be reducible to a 
potential problem, the solution of which is known once we can determine 
the conformal mapping which transforms the exterior of the cross-section 
into the exterior of its ‘equivalent circle’, the modulus of the transformation 
being unity at infinity. ‘This is the problem in electrostatics of determining 
the potential due to an infinite cylinder with a charge distribution of one 
unit of charge per unit length of the cylinder. A(x) for a given cylinder 
is then determined from the momentum integral equation, Glauert & 
Lighthill’s solution for a circular cylinder being a particular case. 

It is shown that, only for flows down cylinders whose cross-sections 
contain a re-entrant angle of 90° or less, does this solution break down 
qualitatively at the leading edge by not giving initially the expected 
square-root law of boundary-layer growth with distance from the leading 
edge. Further for a large class of cross-sectional shapes, including circles, 
ellipses, convex polygons, and partly-concave polygons having re-entrant 
angles greater than 57/8, this method gives quantitative results at the leading 
edge for bulk properties, such as skin friction force per unit length, and 
displacement area, which are quite close to the known Blasius value, 
and much better than values obtained by the Rayleigh methed. On the 
other hand, there is of course no attempt to predict velocity distributions 
in this or any other Pohlhausen method. 

In addition, the results are expected to be good, like those given by the 
Rayleigh method, far downstream, because the assumption of a profile 
for which the viscous retarding force is zero becomes closer and closer to 


2Q2 
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the truth. ‘Therefore, it is not unreasonable to suppose that the complete 
distribution of these bulk properties as functions of the distance from the 
leading edge is given satisfactorily. 

Solutions in closed form are given tor the flat plate of finite width and 
the elliptic cylinder. The skin frictional force and the displacement area 
for the flat plate are given with errors of 3% and 10% respectively at the 
leading edge, and for the circular cylinder errors of 13% and }% in the 
corresponding quantities. For the general elliptic cylinder the corre- 
sponding quantities are given with errors which lie between these two 
sets of values. (Rayleigh’s method gives an error of 70% and 35% in the 
skin frictional force and displacement area respectively.) It is expected 
that the errors gradually become less, further downstream. 

The solution of the problem for the general cylinder is given in the form 
of a double integral which is expanded in series both at the leading edge and 
far downstream. ‘These solutions can, in general, be easily joined giving a 
solution valid along the whole length of the cylinder. The bulk properties 
of the flow are seen to depend on the cross-sectional perimeter at the leading 
edge, and on the ‘ equivalent radius’ far downstream, the detailed shape of 
the cross-section being important in the intermediate region. 

Results for the finite flat plate and the square cross-section are given 
in tabular and graphical form. In particular, the graph of momentum 
defect area enables easy readings of the total frictional drag of the cylinders 
of any length to be made. Some consideration is also given below to several 


other polvgonal shapes. 


(SENERAL THEORY 


If w is the velocity in the x direction parallel to the generators of the 


: ' , ; 1 
cylinder whose boundary is denoted by I’, and z, zw are the velocities parallel 

to Cartesian axes y, 2 fixed in the plane of tne cvlinder cross-section, thet 
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W here 
—- ul uy dvd (5 


is the momentum detect area, % being the whole area outside the perimeter, 


F nd aa (6) 


(where dr, dn are elements of the perimeter |’ and of its outward normal) 


and 


is the frictional force per unit length of the cylinder. Equation (4) may be 


obtained either by direct physical argument or by integrating (1a) over the 


whole area outside | and using (1b) to eliminate v and w. 
The viscous retarding force given on the right-hand side of equation (1 a) 
hand ‘ 1 +4 r\ 4 ; 
vanishes for the profile (3) 1 
2% ; 
— 1}, (/) 
C2” 





The particular solution of (7) with dé = 0 on V which will be used ts simph 
the electrostatic potential of a cylinder at zero potential and unit charge 
per unit length, so that 


Oo 
| —dr = 2. (8) 
y on 


The join between the protile chosen for the boundary layer and the constant 
value U, which wv must take outside the boundary layer, is effected by 


writing A(x) = U//x(x), where x(x) is non-dimensional, and putting 


u= Ud(y,2)/a(x) for 0 --d <4, 
u = l tor db = ot, (YQ) 


Then by (6) and (8) 
F = 2np A(x) ap l//a(x). (10) 


When ¢ has been determined, the variation of x with w follows from (4), 
which, with (8) and (10), gives 


d {| b(a -p) ; dav 





san ee dydz re SER aC (11 ) 

dx ! . Ua(x) 
Once the integral in (11) has been determined as a function of x, It Is easy 
to obtain the relation between « and x by one further simple tntegration. 


\ seemingly objectionable feature of this method 1s the discontinuity 
in du/en at the edge of the boundary layer (given by 4 = x). However, 
Glauert & Lighthill have shown that this in itself does not lead to serious 
errors in a Pohlhausen treatment. Moreover, as vx/(Uc) increases, the 
discontinuity (which is proportional to 1/«) decreases and we would expect 
the errors arising from it to decrease also. 

It now remains to determine the integral on the left-hand side of (11) 


as a function of « for a given boundary shape |’. ‘Che solution to the potential 


problem given by (7) and (8) has been found tor a number of boundary 
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shapes. In particular for a circular cylinder of radius ¢ the solution takes 
the form 


o= log(r/c), (12) 
and if ¢ = re’’ then this can be written as 
cel 10 _ ce’, say. (13) 


The solution of (7) and (8) for a given boundary I then reduces to finding 
the conformal mapping 
w=f(¢), with f(f)~ ¢ for large @, (14) 
which maps the exterior of the boundary I in the w (= 2+7y) plane on to 
the exterior of the equivalent circle of radius c in the ¢-plane. The solution 
of 
w = f(ce*) (15) 
then gives ¢ (the real part of s) as a function of y and 3. 
If we now change variables in (11) from the physical (y, z)-plane into 
the (6, ¢4)-plane, noting that 


dyad dw/ds * ddd, (16) 
the equation becomes 
d ;~d(a—d) r= | dw |? 2var a 
a | eee teen. GAP se ee, (17) 
a5, Jo | ds | Ua 


‘The mathematical problem is thus reduced to the calculation of the inner 
integral in (17). 
A quantity of special interest is the displacement area 
ce 
\=—]|| (U-u)dydz (18) 

UsJs : 
which represents the amount by which the fluid outside the boundary 
layer is displaced owing to the reduced flow inside. On the present theory 
(18) takes the form 


19 





tly —d Ww 
il Li mi ddd, (19) 
which is a function of x. Once (17) has been solved numerically to give 
the variation of x with (vx)/(Uc?), equation (19) enables us to determine 
the variation of A with (vx)/(Uc?). 

To check the accuracy of the method for a given cylinder we expand the 
solution at the leading edge in terms of the parameter (vx)/(U/?) and compare 
it with the expected Blasius layer, where / 1s the perimeter divided by 27, 
given in the present method by 


neg dé. (20) 


3. FLAT PLATE OF FINITE WIDTH 


‘Vhe well-known mapping which transforms the exterior of the finite 
Hat plate of width 4c in the w-plane on to the exterior of its equivalent circle 
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of radius c in the ¢-plane is 


w= C+c/C, (21) 
which can be transformed by using (13) into the required form 


w = 2ccoshs. (22) 





The form (22) shows that this method assumes that lines of constant 
velocity are confocal ellipses when at the leading edge they should be flat 
cylinders. ‘This is a source of error at the leading edge, although the effects 
may well cancel out more or less on integrating round the perimeter. 

On substituting the value of dw/ds* given by (22) it is easily deduced 
that, for a flat plate of width 4c, 

“21 


dw/ds|* d@ = 2mc(e** + e-**), (23) 





which when substituted in (17) reduces the momentum integral equation to 





f (PEE) (cab 4 em) dg = - 


dxjJo Ucta 


The integral in (24) is 





(1/«?)(« cosh «— sinh «)cosh a, 





whence (24) can be solved giving 
Vx 1 . ee =* cosht—1 . 
—3 = 5) cosh2a+3—-sinh2a+ | ————dt}. (25) 

ve «2 a t 
‘The integral in (25) is a well-known exponential integral which has been 
copiously tabulated so that the variation of « with (vx)/(Uc?), and hence 
the skin frictional force per unit length # = 2zpuU/x, is easy to determine 





~U 


numerically. 

For practical purposes the results are given in tabular and graphical 
form. ‘The only means of finding if the solution is accurate up to the leading 
edge is to expand the solution in powers of {(vx)/(U/’)}'* and compare the 
first term in the expansion for skin frictional force and displacement area 
with the expected exact values of Blasius. For the flat plate of finite width 4c, 

l = 4c/z. (26) 

The right-hand side of (25) may be expanded as 

; ({T],;2 1 ~y2 344) a6 
(vx)/( Ue ) re T 39% O(«°) 


which is easily inverted to give, using (26), 


F I { vx ae vx \l2 vx \3/2 vx \5l2 
——$—$<$—$—=<— — — ()- 0, ———- } ° () ace } 8.9 : rn O sud ; 
Fy ~~ ONT p) “a ( B) +! ( F) (air) 

(27) 
‘The exact coefficient of {(vx)/(U/)}-'* in (27) is given by Blasius as 0-332 


showing that even at the leading edge the method of the present paper 
applied to a flat plate gives the frictional force with an error of 3%, although 





some cancelling of errors doubtlessly occurs. Since it is reasonable to 
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assume that the method becomes more accurate tar downstream we see 
that the distribution of skin frictional force is given remarkably well along 
the whole length of the flat plate. In contrast the Rayleigh method gives 
the coefficient 0-564 which overestimates the skin frictional force by 70%, 


at the leading edge. 
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‘l’o determine the displacement area \ we substitute (23) in (19) to give 


A *y—h fet yr, 1 cosh 2x — 1 04 
= | — (Pte) db = 5 ——_—_. (28) 
27c? | , a 2 x. 


Kquation (28) together with (24) enables us to find numerically the 
variation of A with x. ‘lhe results are given in graphical and tabular form 
in table 1 and figure 1. 

For conditions at the leading edge we expand the right-hand side of (25) 


as 

A/(2mc?) = a + 403+ O(a°). (29) 
When we substitute for ~ in (28) as a power series in {(vx)/(UP);\'? we 
obtain 


A vx \lie ye \ 2 ( ye \ i= 
a -Q24 — ). _ 1 ( —. A 30) 
ie 124 ( Tr) : 18 (Tis) (zn) 3 


Blasius gives 1:72 as the exact coefhcient of {(vx)/(U/*)}'*, showing that 
at the leading edge the Pohlhausen method of the present paper gives an 
error of 11% in the displacement area. ‘his can be compared with an 
error of 30% in the Rayleigh method which gives the coefhcient as 1-13. 

In addition to A, the momentum defect area A® defined as in (5) will 
be considered. ‘lhe integral in (5) has already been evaluated to give 


© /(2mc?) = «-*(a cosh « — sinh «)cosh x, (31) 


whence /(27c”) may be plotted against (vx)/(Uc?). 

No comparison with other solutions need be made as © 1s simply 
proportional to the integral of the skin friction / with respect to x. ‘he 
graph of © as a function of x gives a convenient way of determining the 
drag pU?@ on a plate of length x. 


4. ELLIPTIC CYLINDER 
‘The mapping 
w= 6+cRIL (32) 
transtorms the exterior of the ellipse with axes c(1 +k) and c(1—R) in the 
w-plane on to the exterior of its equivalent circle of radius ¢ in the ¢-plane. 
We may assume () k | since this range ot k gives all possible values 
of the ratio of major to minor axes. 
Equation (32) together with the transformation (13) gives the solution 
of (7) and (8) in the form 
w = c(e’+ke-), (33) 
from which it is easily deduced 


“Qn 


dw/ds\|2 dd = 2mc*(e* + k?e-**), (34) 


and 
| = (2c/7)(1+R)E[2k 471 +R) '] (35) 


where F is the complete elliptic integral of the second kind. 
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Substituting (34) in ( 


17 
al pce P) 2 26 + ke) dd = 


), the momentum integral equation becomes 





(36) 


Waa’ 


which is easily solved to give 











For conditions at the leading edge we expand the right-hand side of (37) 
to give 

(vx)/(Uc?) = 3(1 + k?)o? +3(1 — R*)o3 + 2(1 + R2)a4 + O(a), (38) 
which is we inverted to give fails with (35), for the skin frictional 
force F, 




















_—. _ 0-455 (1 +k?)"? / vx \-12 Loner 1—k 
Im (I+R)E aR) TER 
- (1+R)E / vx \!2 vx \3!2 
10-995 Sr Loo) . 9 
> (Rye um) ur) ° ©” 
The coefficient of {(vx)/(U/*)}** is a monotonically increasing function 


of k and varies bias een 0-289 for a circular cylinder (k = 0) and 0-321 for 
a flat plate of finite width (k= 1). Since for any cylinder cross-section 
we expect the error in this method to decrease far downstream it is not 
unreasonable to suppose that for the case of a general elliptic cylinder the 
frictional force is given everywhere with an error less than 13%. 
We can determine the displacement area A in a similar manner to that 
given for the flat plate and find 
p2X __ ) 
sas = at Fx 1+A7[1-—3(1-e ‘Ip. (40) 





‘The right-hand side of (40) can be expanded to give 
A/(2ac?) = $(1+k?)a + 4(1 — R?)a? + O(a), (41) 


which, with the solution of (38), gives 


A _. (14+ Rk2)12 / vx \12 vx \3/2 
772 l ay ale ) = mon O a . (42) 
27/2 1+k)E \UP UP? 
The coefhicient of {(vx)/(U/*)}'? is a monotonically increasing function of k 
and varies between 1:73 for k = 0 and 1-92 for k = 1, showing that it is not 





unreasonable to assume that the displacement area A is everywhere given 
with an error less than 10°%. 


In addition the momentum displacement area © is given as 
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5. POLYGON CROSS-SECTIONS 
We now consider the application of the present method to flow down 
cylinders whose cross-sections are n-sided, closed polygons with exterior 
angles 7a, (r = 1,...2). Suppose that we can determine the a’s of the 
Schwarz—Christoffel transformation 


’ dw/dt = A(t—a,)...(t—a,)™/(1+#), (44) 
with 
Se. = Z, > «,/(a? 4 I) = I, y a, a,/(a-+1) = 0, (45) 
r=1 r=! r=] 


which transforms the exterior of the closed polygon in the w-plane on to 
the lower half of the ¢t-plane, where A is a constant and conditions (45) 
ensure the polygon is closed and that the solution of (44) is single-valued. 
Then introducing the bi-linear transformation 
t =1(€—c)/(f+c), (46) 
which maps the lower half ¢-plane on to the exterior of the circle of radius c 
in the ¢-plane and which with (13) can be written 
t = tan(zs/2), (47) 


we see that the solution of (7) and (8) reduces to solving 








dw z 
— = BB] [sin 3(6+ B,— 41), 48) 
ds = BL [sin 2(0 + B,— gi) (43) 
if the constant B is suitably chosen, where 
a, = tan(38,.) es ee) (49) 
The conditions (45) become 
da, = 2, d«,cosf, = da,sinB, = 0. (50) 
r=1 r=1 r=1 
From (48) it is easily deduced 
dw uf. cos(6+B,)) 
—| = B' cosh <1—- ———*> . 51 
ds eT] | cosh¢ | iad 
Since we seek the solution of (55) with 
dw/ds| ~ ce® for large 4, (52) 


we see that B’ = 2c. 
For conditions far downstream we expand 
dw/ds|? dé (53) 
“0 
in (17) and (19) in a series valid for large @ which can always be written, 
by (50) and (51), as 


on 


| dw/ds|? dO = 2nck(e?+ > A, e-”*). (54) 
0 m=1 


If we use the form (54) the solution of the momentum integral equation 
giving the variation of «, and hence of frictional force, with (vx)/(Uc?) is 


x 2x __ 1 \ 22 el] “A me 
i = i] e 5 -4(¢ Fo ) r _ at | : 2 A, G(a,m), (55) 
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where 
. lf , deerme pm diet , 
G(a, m) —| e@ "*+3—-4——— | —— dt }. (56) 
m E J0 l 


In addition the variation of \ with large x, and hence by (55) with large 
(vx)/(Uc?) is given in the series form 
co 
A/2mc? =(1/4a)(e?*—2a—1)+ > A,, M(a, m), (57) 
1 
where 
M(a, m) = (1/m?x)(ma + e-'"* — 1). (58) 
lor conditions at the leading edge we expand (53) in a series valid for small 4 
and solve equations (17) and (19) in series valid for small (vx)/(Uc?). In 
particular the expected square-root behaviour in the boundary layer at the 
nose occurs if the integral (53) converges for d = 0. 


Regular polygon 
For an n-sided regular polygon it is easily shown (Bickley 1929) that 
3. = mrin ie wee 8 (59) 
whence 
dw/ds| = 2’c\sin(s/A);’, (60) 
where 7A = 27/n is the common exterior angle of the regular polygon. 
From (60) it is easily deduced that 
(27 |dw 
to 
Ja Was 
The integral on the right-hand side of (61) is a hypergeometric function 
which can be expanded in series valid for large and small ¢. 
For conditions far downstream we expand (61) in a series valid for 





“dO = 2'*1¢2 cosh(2¢/A) (1 — sech(24/A)cos y)* dy. (61) 


large ¢, namely 
-2n /? 


d 
| | dO = 2mce*” F,(—A, —A; 1; e~*"), (62) 


which is in the form given by (54) and enables us together with (55) and (57) 
to obtain numerically the variation of frictional force and displacement area 
with large distance x downstream. 

At the leading edge we expand (61) in the form 


27 dw)? | ~ (2A!) 
red ao 2xc7e* Es: F(- A, —A; -2i: l—e aol) 4 


(A!)? j 2) Y/ dl F 
TAD “gyi tan 2A e404) +24 (1 +A,1+A; 242A; 1-e* | (63) 


except for the case of a square cylinder (A = }) when the integral on the 
right-hand side of (61) is an elliptic integral. ‘The expansion of (61) given 
by (63) enables us to determine a series for the left-hand side of (61) valid 
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for small d which when inserted tn (17) and (19) readily gives for non-zero A 


F | C(A/2)! 72 vx \-!2 ye \/ 
Bs hen (2, 178 2f 2X me} Paci 
f° ee | U(2A) 1] (zi) b(t) 


vx \ 12 vx \! 
T Ro } F+TIO T ( > ’ 
ks(a(4 p) (z 7) (64) 
and 


A eal BIT Pavel vey vx \! 
nl? | Al VL | (20) ( T73) t Dt) 
vx \32 | ee dil P 
t DIN Tp) + Ol ap) ‘ (05) 


l= cA$/(A/2)! (66) 


where for a regular polygon 
and 


e=4} A>} (67) 


, 


y 


‘The coefficients of {(vx)/(U/)}"? in (64) and of {(vx)/(UP)}? in (65) vary 
between the limits given for a flat plate (A = 1) and a circular cylinder 
(A = 0), whence it is not unreasonable to suppose the skin frictional force 
per unit length and the displacement area are given everywhere with errors 
less than 13° and 10%, respectively. 

An indication is given in (64) and (65) of the effect of flow outside an 
angle m(1—A) on skin frictional force and displacement area at the leading 
edge. For small values of A the coefficients of the fourth terms in (64) 
and (65) are so large compared with the second and third coefficients that 
the contributions of the second and third terms are well below the level of 
the fourth. For values of A~ } the second and third terms effectively 
cancel in both (63) and (64), while for values of A ~ 1 the third terms give 
the effective first approximation. 

‘The results given so far for skin frictional force agree qualitatively with 
a result given by Batchelor, namely: ‘he friction on cylinders (with 
inward curvature everywhere) with corners is always less than the friction 
on cylinders with smooth boundaries; to turn a corner sharply is more 
economical on drag than to do it gently with the same perimeter”. ‘lhe 
results do however disagree quantitatively with Batchelor by not predicting 
a finite effect for corners at the nose. 

In addition a similar result holds for displacement area which can be 
stated: ‘“‘'The displacement area on cylinders (with inward curvature 
everywhere) with corners is always greater than the friction on cylinders 
with smooth boundaries; to turn a corner slowly is more economical on 
displacement area than to do it quickly with the same perimeter”. 

‘lhe values of F'(27U) for a cylinder of square cross-section are plotted 
in figure 1. For (vx)/(Uc?) | they are indistinguishable from those for 
a flat plate or circular cylinder of the same equivalent radius ¢. Glauert & 
Lighthill’s recommended curve, obtained by adding 9% on to the Pohlhausen 


values, is probably a close approximation for all shapes in this range. 
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n-pronged cross 
To investigate the application of the method of this paper to flows down 
cylinders whose cross-sections contain re-entrant angles we consider the 
case of the cylinder with a cross-section formed by » equal, equally spaced 
prongs, whence the common re-entrant angle is 27/n = 7A. The solution 
of (7) and (8) is such that 
dw/ds* = 2’c sin(s/A)]*~ | cos(s/A)|, (68) 
whence J = 2***¢/(aA) (69) 


and 
| |dw/ds!* dO = 2mcte**[ F(—A, —A; 1; e-*” \+ 
+ 4(1 —Aje~*4!* F(2—A, 1—A; 2e-**)]. (70) 


The right-hand side of (70) can be written 
ont > Aa, (71) 


r=0 
where A, = (A— 1)? —(A—r4 1)?(A— 2r)?/(r!)2, (72) 
which is in the form of (54) and enables us to find the distribution of skin 
frictional force and displacement area far downstream. 

‘lo obtain a check on the accuracy of the method at the leading edge we 
continue (71) analytically for small ¢ and obtain 


2 2a)! 
| ta] = 2ne| TOF + Be? 0) |. (73) 


‘The right-hand side of (73) shows that only for re-entrant angles greater 
than 90° (A < 3) does the method of this paper yield results which are 
qualitatively correct at the leading edge since they give the expected 
square-root growth of boundary layer with distance. For re-entrant angles 
less or equal to 90° the surfaces of constant velocity assumed by this method 
at the leading edge are so different from those in reality that the momentum 
defect is seriously overestimated, whence the skin frictional force is over- 
estimated. However, for re-entrant angles lying between z (flat plate) and 
57/8 the skin frictional force is overestimated by between 3°% and 21% 
while the displacement area is overestimated by between 11°, and 40%. 
‘These results are better than those given by the Rayleigh method. 
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Cavitating flow about a wedge at incidence 
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SUMMARY 
A mathematical model is constructed for cavitating flow past 
a wedge with sides of equal length but with its axis of symmetry 
placed at an angle to the incident stream. ‘The model involves 
a subsidiary cavity with a re-entrant jet at the vertex. Only the 
case of zero cavitation number is considered. ‘The flow field 
is worked out in some detail for small angles of incidence, and the 
lift, drag and moment coefficients are calculated as far as first-order 
terms in the angle of incidence. It is shown that the effect of the 
rate of loss of momentum in the re-entrant jet on these force 
coefficients is negligible to this order. 
Experimentally, it is shown that the secondary cavity does 
exist under suitable conditions, and the force coefficients 
obtained agree with the theory. 


1. INTRODUCTION 

The problem of two-dimensional steady cavitating flow with zero 
cavitation number past a symmetrical wedge in an infinite stream, neglecting 
external forces, has in the past been solved exactly, using complex variable 
techniques, for the particular case of a non-yawed wedge (see Lamb 1932, 
or Milne-Thomson 1949). 

Milne-Thomson has also derived a solution involving a yawed wedge. 
He found, however, that for the flow pattern considered it is not possible 
to specify arbitrarily the ratio of the lengths of the two wedge faces. ‘This 
ratio is determined as a function of the wedge and yaw angles, and is unity 
only for zero yaw; so that, in particular, Milne-homson’s method fails 
to solve the problem of a symmetric wedge at non-zero yaw. 

To be able to specify the ratio of the lengths of the wedge faces 
independently of the wedge angle and yaw angle, it is necessary to relax 
the conditions determining Milne-‘Thomson’s problem, and it is clear 
that the non-essential condition imposed is that which requires the vertex 
of the wedge to be a stagnation point. ‘This is borne out by considering the 
special case of a yawed semi-infinite flat plate, whereon the stagnation 
point is not in general in a symmetric position (see Milne-'Thomson 1949). 

It is by discarding this condition that the following solution is found, 
involving a flow pattern which is of a sufficiently novel type to appear to 
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inerit description in some detail. ‘The novel feature of the flow —a subsidiary 
cavity in the lee of the wedge vertex--is shown experimentally to exist 


under certain conditions. 


2. ‘THE PHYSICAL MODEL OF THE FLOW 

The method of solution used will be that of the complex variable, and 
for this we need to assume potential flow and neglect viscosity. 

It is convenient to adopt the convention that the angle of incidence is 
positive, so that the ‘lower’ face of the wedge is inclined to the stream at 
a greater angle than the ‘upper’ face. By analogy with the behaviour in 
the case of the flat plate, we should expect the stagnation point to move to 
the lower face of the wedge. 


2b 





a 2c 
ng 


Figure 2. 

Vhe dividing streamline through this stagnation point will, on the 
downstream side, follow the wedge face until it reaches the rear of the 
wedge and will then become the lower free-streamline of the main or ‘ rear’ 
cavity. On the upstream side of the stagnation point the dividing streamline 
will reach the vertex of the wedge and would then, in the absence of a 
stagnation point, need to attain an infinite fluid velocity to turn round the 
sharp corner of the wedge. However, owing to the inability of real fluids 
to sustain pressures less than vapour pressure, a rise in velocity above a 
certain finite value, in this case taken to be the main stream velocity, is 
impossible. We shall assume in our model that a cavity will form, bounded 
by a free-streamline which is the smooth continuation of the dividing 
streamline. Besides the rear cavity, therefore, we have a ‘ bubble’ cavity 
formed on the upper face of the wedge; this certainly exists under suitable 
conditions, as is shown by the photograph (figure 1, plate 1) taken in the 
\.R.D.E. Cavitation ‘Tunnel. 

The free-streamline springing from the vertex must either rejoin the 
upper face of the wedge or else continue to infinity. ‘lhe first case might 
conceivably occur as illustrated in figure 2. If, however, the re-attachment 
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Figure 1. 30) semi-angle wedge. 


Angle of incidence 10 , cavitation numbe1 0-003 
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were effected tangentially, as in figure 2(a@), the free-streamline would 
contain a point of inflection, which is impossible because the pressure 
gradient, and therefore the streamline curvature, must always be directed 
towards the cavity. If the re-attachment were other than tangential, as 
in figure 2 (b), it would involve a stagnation point, which again is impossible 
on a free-streamline. 

If, on the other hand, the free-streamline continues to infinity down- 
stream, so that the upper face of the wedge is included in the cavity, 
as in figure 2(c), then the flow is simply that past a yawed flat plate. 
This flow pattern will indeed occur under suitable conditions, which 
may easily be defined from the exact general solution of this type of flow 
(see Milne-’Thomson 1949). 

When the dimensions and orientation of the wedge preclude the 
possibility of ‘yawed flat plate’ flow, we are left with the possibility that 
the flow continues to infinity in some other direction. ‘This apparently 
unreal configuration, involving a double-sheeted flow pattern, is similar 
to the familiar re-entrant jet models occurring in the solutions obtained by 
Kreisel (1946) and Shiffman (1949) for non-yawed wedges at non-zero 
cavitation number. 

While Gilbarg & Rock (1946) have shown by numerical methods that 
drag forces are remarkably insensitive to the assumptions made as to the 
flow pattern at the rear of the cavity, at least when this region of the flow is 
far removed from the wetted nose of the body, we cannot quote this as 
sufficient justification for the model now proposed, since the re-entrant 
jet is here postulated as forming on the face of the wedge itself. ‘The 
justification will rest partly on the degree of agreement achieved with 
experimental results, and partly on the mathematical demonstration that 
the rate of loss of momentum in the re-entrant jet is of smaller order of 
magnitude than the drag and lift forces on the wedge, at least for small 
angles of yaw. 

We now consider the assumed mathematical model in detail. Referring 
to figure 3, we have a symmetrical wedge DCD’ of semi-angle § in an 
infinite stream which makes an angle « with the axis of symmetry of the 
wedge. As is customary, the uniform stream is considered to be due to 
a source system at A, the point at infinity (A,, Ag, etc. all refer to the unique 
point at infinity in this plane and are labelled differently for the purpose of 
identification of individual streamlines). 

The streamline A,B divides at the stagnation point B on the lower 
wedge face CD’. A, BD’ becomes the lower bounding free-streamline D’ A, 
of the main cavity and A, BC breaks away smoothly from the vertex C to 
form the bounding free-streamline COF of the subsidiary cavity. ‘The 
free-streamline COE and part of the flow to its left turn almost completely 
round and flow to infinity as a separate jet, bounded on one side by the 
extension of the wedge face DC and on the other side by the free-streamline 
COE. ‘The point at infinity on this jet is denoted by E and is quite distinct 
from A, the point at infinity in the main stream. Since the streamlines 
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in the jet cut across the streamlines of the main flow pattern, it is necessary 
to visualize the flow as occurring on a double-sheeted Riemann surface, 
one sheet of which carries the main flow pattern, while the jet crosses on to 
the second sheet. ‘The cut dividing one sheet from the other is somewhat. 
arbitrary, but it can be said to have the form GOPA, shown in figure 3. 
Since part of the flow to the left of A, B flows upstream and the rest flows 
downstream, forming the upper wall of the main cavity, there must be a 
stagnation point H on the upper face CD of the wedge. ‘The dividing 
streamline A; H will, on the downstream side, become the upper bounding 
free-streamline DA, of the main cavity. 











PHYSICAL PLANE] 2“ RIEMANN PLANE 
WEDGE FACES —_— — 
|ISTREAM LINES rs — «x 
FREE STEAMLINES) ——-——-——- | —-—-—-— 


IRIEMANN CUT __| “ _| 























Figure 3. Assumed model of the flow. 


‘Thus, of the fluid flowing from A, that lying to the right of A, B returns 
to A, and that to the left of A, B and to the right of A; H flows to E and 
is ‘lost’, while that to the left of 4; H returns to A. 


3. "THE SOLUTION OF THE PROBLEM 

We use the standard technique of obtaining conformal transformations 
between suitable complex variables to determine the flow pattern and all 
physical quantities. 

lor the physical plane as represented in figure 3, we use the variable 
z=x+1y, where x and y are rectangular Cartesian coordinates parallel 
and perpendicular respectively to the plane of symmetry of the wedge, the 
direction of increasing x being downstream. Next, we define the complex 
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velocity potential w= 46+, where ¢ is the velocity potential and % the 
stream function. ‘The velocity vector (u,z) is given by 
od Cus od ous 
u — = —. Z — —, 
Ox y oy Ox 
By a suitable choice of scale in the w-plane, in respect to the as yet 
undetermined length scale, we arrange that 
u=, +77, = 1, 
where (u,,,U,) is the velocity vector of the flow at infinity. Then 
dw/dz = u—iv = ge, 
where g? =u?+v" and @=tanv/u. 

Ihe w-plane is illustrated in figure 4+. All the streamlines, since they are 
by definition lines of constant 4, become lines parallel to the ¢-axis. ‘The 
pair of streamlines HDA, HE and the pair BCE, BD'A should each actually 
be represented by a single straight line but are shown separated for 
convenience. 

A H D Ag 


5 ——-—-— — —-—-—e 
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Figure +. ‘The complex velocity potential. 


(‘he next variable is the logarithmic hodograph variable Q defined by 
Q) log(dw dz) log g—10. 

This variable is plotted in figures 5(a) and (). ‘The free-streamlines D’ AD 
and COE are lines of constant speed q, this constant speed being equal to 
that at infinity so that, by our choice of scale, we have g = 1 (i.e. logg = 0) 
on the free-streamlines. ‘Thus in the Q-plane they become portions of the 
imaginary axis. Since the wedge faces CB, BD’ and DH, HE are lines of 
constant #, they become lines parallel to the real Q-axis with a discontinuity 
of amount z in the imaginary part of Q at H and B, where Re(Q) = — oo. 

It is found on tracing the other streamlines that the Q-plane is doubly- 
covered. ‘l’o overcome this we make a cut in the 2-plane from the point R, 
where d{2/dz vanishes, to the point at infinity in the Q-plane in a negative 
direction, and then consider the {2-plane as a Riemann surface with two 
sheets cross-joined along this line. 

‘Then the streamlines are as shown in figures 5(a) and (6), where, for 
convenience, we have separated the two Riemann sheets. ‘lhe streamlines 
themselves show how the four half-planes are joined. ‘lhe images of all 
streamlines between the dividing streamlines A, B and A. H ot figure 
lie between AX, and AX, in figure 5(a); after crossing the cut they remain 
on the second sheet, eventually termfnating at F in figure 5 (bd). 
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In order to obtain the requisite conformal transformations between 
these variables, we define an auxiliary variable + by requiring that the 
streamlines DAD’, D'BC, COE, EHD become the sides of the rectangle 
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Figure 5. ‘The logarithmic hodograph. 
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Figure 6. (a) 7-plane, (0) t-plane. 
formed by the four points r= +K, +K+7K’, the orientation of the 


lines being as shown in figure 6(a). We also suppose the points H, B, A 


to be represented by 


RK’ >a,5 20 


7, = 1K' +a, -K <a<kK. 
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In this plane the streamlines are as shown in figure 6(a), the behaviour 
near A being precisely as in figure 5(@) apart from a change of orientation. 

We now find by inspection the conformal transformations which 
transform the lines DAD’, D'BC, CE, EHD and AB, AH into their images 
in the w- and Q-planes. From this we shall find a parametric relation 
between w and dw/dz, and hence the problem is solved. 


The Q-r transformation 

‘The functions to be used in this and the following transformation, as 
has been foreshadowed by choice of the 7-plane, are the Jacobian elliptic 
functions, for the definitions and properties of which reference may be 
made to Neville (1951), Jahnke & Emde (1948), Milne-homson (1950), 
and Whittaker & Watson (1940). Although Neville’s notation is certainly 
the more elegant, we have used that of Whittaker & Watson as being the 
more familiar. 

The function dQ/dr must have a simple pole at ry and 7, and must also 
be purely imaginary all along the boundary rectangle in the r-plane. It 
must also be finite and non-zero at 7,. ‘The transformation 


dQ) T,snt T.snt 





pidicie —- -alV. (1) 
dr snt—Sntg sntT—SsnTp, 
where N is real and 7, 7, are purely imaginary, satisfies all these require- 
ments since sn7 is real all along the rectangle in the r-plane. ‘The function 
dQ/dr can be shown to have only one zero inside the rectangle and this 
must obviously be at tp, where R is the branch point illustrated in figure 5. 
To ensure that the transformation gives the correct scales as well as 
directions we have to impose further conditions. ‘The first two are that Q 
must have a discontinuity of i7 at ty and 7,, and these are easily seen to 
lead to 
T', = cntgdsrg, and 7, = cntzdsrTp. 
The other conditions are 
QQ» ip, Q) 12, Qa -1(7—), QQ), i(z — B), Q, - 10, 
and these must be applied to the integrated form of equation (1). Integrating, 
we get 
Q). = I(r, 77 —7K’) + Il (7, 73-1 K’) 
bhogel ee ee ie ie 
oe Cel TH COT cd Tp cdr , 
where ‘ck? snyenydnvysn? u 
I(x, y) Biles aes oem ee 
I R°sn* ysn°u 


as defined by Whittaker & Watson (1940). ‘The constant of integration 


in (2) is zero due to the skew symmetry of the conditions at C and E. 
Che first four conditions imply that 
NK = (x7—8)+i[II(K, 7, —1K’) + I(K,7,3—1K’')}, 


SK’ 
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The fifth condition gives 
-ta = II(a,7q —1K") + l(a, 7, —iK’) + tatZ’(h) + 2'(b)} - 
( ksc’b ksc’h) 
ta Se eee 
dca dea } 


if 


; (  (k’)? snasn’hen‘h) 
retn + — - r 3} 
cnadnadn’h | 





; if (k’ snasn’bcn’b\ 
+2 tan~< 





7 "5 4 
cnadnadn’b f{ *) 
where the prime indicates that the co-modulus k’ replaces the modulus k, 


and Z(x) is the Jacobian zeta function. 


The w-z transformation 

This transformation is more easily found by considering a further 
auxiliary variable ¢ defined by t=snzr. Since snz is real all along the 
rectangle in the 7-plane, the ¢-plane is as in figure 6(d). 

Then, using the Schwartz—Christoftel theorem, we see that we require 
the following conditions for dw/dt: 

(i) It should have a simple zero at ty, tg. 

(ii) It should have a simple pole at ¢, and a triple pole at ty. 
Thus we get 

dw P(t—ty)(t—tz) 


dt = (t—t,)(t—ty)’ 





or, in terms of 7, 


dw Pcnz7dn7(snt—snty)(sn7—sn7z) (5) 


dr (snz—1)(snz7—sn7,)° 





since sit, = snK = |. 
The constant P determines the length scale in the problem, and it is in 
relation to this that we choose the scale in the w-plane so that q, = 1. 
Equation (5) may be integrated to give w as a function of 7. Also, since 
Q2 = log(dw/dz), equations (2) and (5) together yield an expression for 
dz/dz which can again be integrated to give z as a function of 7. ‘Thus we 
can determine w and s parametrically in terms of 7, and hence find the 
flow pattern. From these we can determine all the physical quantities 


required. 


4, DETERMINATION OF PARAMETERS 
‘There are four parameters a, b, h and k in the problem. Equations (3) 
and (4) show that specifying the wedge semi-angle and the angle of yaw will 
provide two equations. ‘lhe other two conditions to be imposed are the 
following : 
(i) The lengths of the sides of the wedge must be in a predetermined 
ratio. In particular we shall consider a symmetrical wedge. 
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(1) In general C will not lie on the continuation of DH. Ensuring that 
it does provides the fourth condition and makes the solution 
determinate. 

Conditions (i) and (ii) can be expressed in the mathematical form 


DD’ = 0+27\/CD'\sin B, (6) 


where DD’ and CD’ are conveniently calculated by integrating dz/dr 
along the sides of the rectangle in the 7-plane. ‘The real and imaginary 
parts of this equation, together with equations (3) and (4), provide the 
four equations necessary for determining the four parameters of the 
problem. 

In general these quantities will have to be determined numerically, but 
in certain cases, and in particular when k is small, the calculations can be 
carried out analytically as far as first-order terms. Since it can be shown 
that small k combined with certain other assumptions corresponds to small 
angles of yaw, a situation of particular interest, we shall consider this case 


in more detail. 


5. ‘THE LIMITING CASE OF SMALL ANGLES OF YAW: THE FLOW PATTERN IN THE 
LARGE 

As the yaw tends to zero, it is reasonable to suppose that the flow pattern 
studied here will tend to the familiar symmetrical flow pattern. In other 
words, the effect of the existence of the secondary cavity on the main cavity 
should tend to zero as the yaw tends to zero. Similarly the secondary cavity 
must become so small that its detailed geometry is unaffected by the main 
cavity. ‘lhus the two free-streamlines must tend to become independent 
of one another, or, interpreting this in the t-plane, K’/K — © as the yaw 
tends to zero. From the theory of Jacobian elliptic functions, this implies 
that k tends to zero. 

Now, writing 28/7 = e, where 0 < e 1, equation (3) becomes 


h+b=eK’, 
and, since h, b > 0, 
h<eK’, b<eK’. (7) 
Writing o = i1K’—7, we transform the two fundamental equations (1) 


and (5) into the form 


dQ cnoy ds oy cnopdsoz 


iN, (8) 





do l—nso,sno 1—nsogsno 


and 
dw Pksn'o,cnodno(1—nsoysno)(1—nsogsno) 


do 
When k -> 0, and provided o lies on or near DAD’, the following relations 
hold: 
sno = sina[1 + O(R?)], cno = coso[1 + O(k?)], dno = 1+ O(R?). 


(9) 





(1+sno)(sno,—sno)? 
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Moreover, it is easily shown from (7) that 
nso, = O(R"), nsog = O(k"), 
cna, dso, = —1+O(k*), cno, ds oy -i + O(k*), 
where l>an=1-—e > 0. 
Equations (8) and (9) therefore reduce to the form 
dQ/do = M—ipsine+ O(R’), (10) 
and dw S cos o 


ne a perenne Ti+ maine +O i, 
do (sino 4 sina! fA sina (/ )] (11) 


where w= —Nsog—nsog = O(k"), 
and M, S, y are constants, y being chosen as the smaller of the two 
quantities 1 or 2n. 
Integrating (10), and applying the conditions at D and D’, we obtain 
>) a aay ‘ : 
Ct = e 1 —tpCOS 0 + O(R’)]. (12) 
Since exp(— 2) = dz/dw, we find that 


dz cosa 





do (sino + sin a) 
Hence 
pate COS Go 
DD'{1 + O(k’)| Pe ‘) snpcorcnsninersnss [e-*** — tue” °| do, (13) 
: »(sino+sina) 
the path of integration being the straight line from — 7/2 to 7/2, indented 
on the positive side by a small semi-circle described about o = —a. 
From the condition (6), expressing the equality of the wedge faces, it 
follows that Re(DD’) = 0. Applying this condition to (13), we find that 
a=yp/t(e) = O(k”), 


the derivation of this formula, and the definition of f(e), being given in 
Appendix |. The function ¢(e) varies monotonically from 4/(4 — 7) = 4-660 
at « = 0 to 5-000 at e 1 (as shown in table 1, $6). 

When this condition is satisfied we have 


DD’ S sin(dez)F(e)[1 + O(R’)], 


where Fe) l+e- se"[h(1 —€ 4) - ub( 1 2- € 4), s(x) being the digamma 
function (:4(x) = (d/dx)V'(x)). In particular, F(0) = 1 and F(1) = 2472/2. 

In order to determine the forces on the wedge, let us consider a circle 
in the o-plane, given by o a+re” (r small). In the g-plane this corre- 
sponds to a large contour ¢ ympletely surrounding the wedge and enclosing 
most of th is shown 1n figure 

hen, equating the forces acting on this contour to the rate f chang: 
ton entum inside ve | 

V,} » | (n. q)qds | onds+R, 

where (X,}) is the force per unit length acting on the wedge, n ts th: 

vector along va rmal to the contour, s measures arc length 
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along the contour, p is the density, p is the pressure, q is the velocity vector 
(u,v), and R is the contribution of the fluid in the re-entrant jet. 

Referring to figure 7, we know that on OSP, q=0, p=p.; and on 
PRQ, q is known, p = p,, + 5p(1—q*). Hence 


(X,Y)=p] (m.q)qds—tp | (1—q*)nds+R, 
“A 


where A is the semicircle o are" 0 n TT. 





Figure 7. 


Since the left-hand side of this equation is independent of r, so also is 
the right-hand side. ‘Thus, on letting r tend to zero, we find 
X = teaSp+ R, + O(a’), 


and } lampS| t+e«(2t—1)—et4 e3] + R, + O(a*), 
where R, and R, are the components of R. 
From (12) we find, on putting o a, that the angle of incidence 
‘ ~ ¥] Ta ) y ; y 2 
is given by z= ale—t' + O(a?). 


Hence the drag force D and the lift force L per unit length are given by 
D XY + O(a?) le*apS + R, + O(2?), 
and 


L= Y—Xxa+O(#?) = sampS|[—t+(2t—1)]+ R,—«R, + O(2?). 


4 
In a similar manner we can find the pitching moment / per unit length, 
measured about the mid-point of the base of the wedge, and reckoned as 


positive for a destabilizing moment. Using the formula 


14 


V/ | p(xdx+ydy)~—p | (vx—uy)udx—vwvdy), 


we find that 
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where R, represents the term due to the re-entrant jet, and the function 
p(e) is explained in Appendix III. In particular we have 
p(0) = 3/(8—27), and p(1) = —3/2. 

The terms R,, R, and R;, representing the effect of the re-entrant jet, 
must be small if our mathematical model is to be physically realistic. In 
the next section we shall prove that they do, in fact, become negligible 
compared with the other terms in D, L and M as « tends to zero. 


6. 'THE FLOW PATTERN IN THE NEIGHBOURHOOD OF THE VERTEX AND ITS 
EFFECT ON THE FORCES ACTING ON THE WEDGE 

‘The analysis of the preceding section was effected without making any 
assumptions regarding the behaviour of 5 and / as k tends to zero, apart 
from the restriction given in equation (7). It is now necessary to examine 
their behaviour in more detail, and it is clear that three possibilities arise : 

(i) 6 and / both tend to infinity; 

(ii) 6 remains finite; or 

(iii) A remains finite. 

The correct choice is gov erned by the tact that we have yet to satisty 
the condition imposed by the imaginary part of equation (6), representing 
the fact that C must be collinear with D and H. (The real part of this 
equation has already been satisfied by the condition x = a t(e) in the previous 
section. ) 

Since the analysis differs in detail for the three different possibilities, 
we shall reproduce here only that corresponding to case (111) above, since, 
as we shall show, this is the case that enables equation (6) to be satisfied. 

In this case it follows from (3) that 6/K’ — e, and hence 

NS Tp ™ ncTz = O(R'*), and dnt, = 1+ O(k?”). 
hen on or near CE it can be shown, from (1), that 
= scot eo O(R*), 
at sIn7— SINT} 
where Z 1s a real constant. 
Integrating, and applying the conditions at C and E, we find 


; cos 3(7 +77) - 1 
f ] a : é [1 + O(k )}- 
sin 3(7— Tz) 
Similarly, trom (5), 
lw 4§SRk!*"cos 7(sin 7 — sin 71) 
teen FE OR. 
dr 2(sin 7 — 1) 
Hence 
dz Skirn 1+sin7 : = : 
+ 14 —— - — |1+cos(r-+ Ty )|{1 . O(R‘)). (15) 
dt a COS T ; 


Now the imaginary part of equation (6), implying the collinearity of C, 


D and H, may be replaced by the following condition: 


fm CEe-™ > = 2, 16) 
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where E is the point at infinity on the free-streamline of the re-entrant jet, 
and g is the width of the re-entrant jet at infinity. 
But g may be readily evaluated by studying the behaviour of equation (15) 
in the neighbourhood of 7 = K = 7/2+ O(k); thus 
g = 3SR''"a[1+sinzg][1 + O(A)]. (17) 


Substituting in (15), our condition becomes 


ot 


m[{cosh A — 1][1 + O(*)] | (1 + cos ¢)(sin t)“"(sin nt 
0 


coshhcostsinnt—sinhhsintcosnt) dt, (18) 
where we have written 7 = 7/2—t. 

In Appendix II this equation is solved explicitly and then, using 
equation (15), the dimensions of the subsidiary bubble are calculated. 
However, as we shall now show, the effect of the value of # on the forces 
is negligible. 

Equation (17) gives the value of g and it is obvious without actual 
evaluation that the values of R, and R, will both be equal to multiples of g, 
the multipliers being of order unity as « tends to zero. However, g itself 


is of order «!+!”, so that we have 


D le2z70S + O(a), (19) 
E lampS[—t+ €(2t—1)]+ O(2?), (20) 
V lap7S” p(e) + O(a”). (21) 


We now define a drag coefficient and lift and pitching moment slope 
coefficients, using the dynamic head as the reference pressure and the 
combined length of the wedge faces c= | DD’ cosecf as the reference 
length. An asterisk will be used to distinguish these coefficients from those 
which would be obtained by taking DD’) itself as the reference length. 
The advantage in taking c as reference length is that all these coefficients 
tend to finite limits as e tends to zero. 

Using an obvious notation, we obtain 
- D €"77 
logz, € F(e) 


oC? E T t+ «(2t—1) 

—# hm | ——; : : ] . (20’) 
eh el +20 L. 3PGa ae («—t) Fe) 

ra Dies - M =e : 

= = lim] ; meld =| = (<) . (21’) 
O@ 1. .6 20 LoP@e ae («—t)F*(e) 


It is interesting to compare these formulae with known limiting values 




















as « tends to zero and as f tends to zero or 7/2, bearing in mind, however, 
that the foregoing analysis is not strictly valid when § equals zero or 7/2 
For the case « equals zero (19’) agrees with the formula given by Perry 
(1952), while (20) and (21) give zero values for L and M as would be expected. 
When « tends to unity, corresponding to a flat plate nearly perpendicular 
to the stream, (20’) agrees with the value given by Milne-‘Vhomson (1949), 
and (21’) with that given by Lamb (1932). 
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When « tends to zero, corresponding to a flat plate nearly parallel to the 
stream, we find that (20’) gives a value z and (21’) a value 37/8. ‘These do 
not agree with the classical values for cavitating flat plates, but this is to be 
expected since even in the limiting case both sides of the wedge are wetted. 
A more reasonable comparison is with a flat plate at incidence in fully 
wetted flow, with circulation chosen using the Joukowsky hypothesis of 
a finite velocity at the trailing edge. In this case Ramsev (1947) gives the 
value z for [0C#/0x], — 9 and 3z/8 for [eC7,/ex], — , so that again agreement 
is secured with classical theory. 

Table 1 gives values of Ch, [0CZ/dx], 9 and [eC},/dx], _y and they 


are plotted in figures 8 and 9 respectively. 


fee | ae: 

, | t(e) Cp | SE | rT | Fal 0 
| | 
| | 





4-660 O 3°: 
1 














0-2 4-721 0-102 I 0-553 
0-4 4-785 0-325 0-629 0-1-7 
0-6 4-853 0-571 | 0-138 =| 0-053 
0-8 4-924 0-770 | 0-628 | 0-110 
| 1-0 | 5-000 0-880 0-880 | 0-092 
Table 1. 
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Figure 8. Modified lift slope coefficient as a function of /. 


7. EXPERIMENTAL METHOD 


1 


\s a result of the foregoing theoretical study of cavitating How about 


a wedge at incidence, experimental data were needed to verify the assumption 














Cavitating flow about a wedge at incidence 629 
made concerning the existence of the subsidiary cavity formed on the upper 
face of the wedge, and also to measure the hydrodynamic forces acting on 
the wedge. ‘These data were obtained in the A.R.D.E. cavitation tunnel 


(see Clayden 1954). 


Bo elz0 “SP 








-02 
Figure 9. Modified moment slope coefficient as a function of f. 


This tunnel is a re-circulating type with a rectangular 9 x7 in. 
free jet working section with fixed vertical side walls and free top and 
bottom surfaces. ‘The working section speed is nominally constant at 
40 ft./sec, and the cavitation number is varied by varying the pressure. 
Very low cavitation numbers (~ 0-01) may be obtained in this manner. 
The basic instrumentation consists of the following: 

(1) a horizontal beam manometer to measure the pressure across the 

contraction nozzle, which is used as a Venturi meter; 

(2) a tilting dead weight manometer to measure the difference between 
the ambient and cavity pressure ; 

a tilting dead weight manometer to measure the ambient pressure 


(3 


— 


in the working section. ‘This instrument is also used to operate a 
servo-mechanism to control the tunnel pressure in the range 
1 atmosphere to 10 mm of mercury ; 
(4) a thermocouple thermometer; and 
(5) a mechanical three-moment balance. 
(‘The first four instruments give continuous recordings.) 

The 30° semi-angle wedge shown in figure 1 (plate 1) was originally 
used to obtain design data for struts and model supports. It completely 
spanned the tunnel, but was attached only to the rear boundary wall of the 
working section. It was used to obtain photographs of the flow pattern, 
in preference to the family of wedges used to obtain force data, because it 
permitted an uninterrupted view. 
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The size of the models used to obtain the force data were as follows: 
tor the wedge semi-angle 8 = 15°, the width of the base of the wedge was 
0-5 in.; for B = 30°, 45°, 60° and 90°, the width of the base was 0-3 in. 
To obtain the force measurements, the wedges were pivoted on the front 
boundary of the free jet and were attached to the three-moment balance 
by a universal joint (figure 10). ‘The three-moment balance was mounted 
on the back of the working section. (A more detailed description of the 
experimental rig was given by Clayden (1957).) A pressure tap was fitted 
in the boundary wall of the working section to measure the cavity pressure 
directly, since the cavity pressure is normally higher than the vapour 
pressure, due to a small quantity of air leaking through from the balance. 


~~ _ NDARY 
iw BACK BOUND 







TRANSPARENT 
FRONT BOUNDARY 
OF FREE JET 


MOUNTING FOR PIVOT 


Figure 10. Test rig. Dimensions in inches. 


For a typical series of tests the cavitation number was held constant 
whilst the balance was rotated from —15° to +15° with the horizontal, 
readings of the three moments being taken at intervals of 3°. Since the 
angle the free jet makes with the model is not precisely known, the position 
of zero incidence was determined as the position at which the measured 
lift force changed sign. 

Although the tunnel is capable of achieving a cavitation number as low 
as 0-01, the lower limit in these tests was determined by the length of the 
cavity, which for the lowest cavitation numbers extended throughout the 


length of the working section. 


8. RESULTS 
Observations and photographs of the flow about a yawed wedge indicated 
that the mathematical model, illustrated in figure 3, closely resembled a 
time average of the actual physical flow (figure 1, plate 1). ‘The subsidiary 
cavity, however, oscillated considerably with a high frequency (a similar 
cyclic process has been observed and studied by Knapp (1956) for flow 
along a cylinder with a cavitating head). As the angle of yaw of the wedge 
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is increased so the subsidiary cavity becomes larger until it finally breaks away 
leaving the leeward face completely unwetted. ‘lhe angle at which this occurs 
has been studied by Clayden (1957) for a family of wedges of varying vertex 
angles. ‘The force measurements given in this paper are for small angles of 
incidence with the leeward face at least partially wetted (figure 3). 














oC, Cy 
p 9 Cp d a na 
15 0:067 0:291 0-00742 0-128 0-157 
0-096 0-300 0-126 0-163 | 
0-119 0-296 0-124 0-159 
0-170 0-345 0-130 0-163 
0-230 0-361 0-141 0-172 
| 
| 
30 0-059 0-517 0-0423 0-0182 | 
| 0-073 0-523 0000305 0:0389 0:0172 | 
| 0-112 0-536 0-0436 0-0192 | 
0-134 0-548 0-0435 0:0198 | 
| 0-183 0-590 0-0438 0-0198 
| 0-206 0-608 0-0431 0-0188 
45 0-081 0-681 0-000342 0-0105 0-00180 
0-116 0-709 0-0097 0-00150 
0-149 0-727 0-0098 0-00212 
|. 0-180 0-755 0-0088 0-:00180 
| 0-260 0-794 0-0097 0-00242 
60 0-116 0-843 0-00032 0-00501 0-00185 
0-122 0-864 0-00523 0-00214 
0-164 0-890 0-00517 0-00248 
0-210 0-885 0-00578 0-00302 
0-268 0-942 0-00488 0:00240 
90 0-104 0-945 0-000170 0-0171 0-00214 
0-130 0-974 0:0174 0-00166 
0-145 0-991 0-0170 -():00148 
0-180 1-025 0-0171 0-00165 
0-220 1-056 0:0168 0-:00231 


























Table 2. 


‘lhe force measurements obtained from the three-moment balance were 
analysed to give Cp, C,, and C,,, the drag, lift and moment coefficients 
respectively. ‘These results were fitted by the method of least squares to 
give the following relations for each combination of vertex angle and 
cavitation number: 

C, = (0C,/0a)a, Cy = (0Cy,/0x)a, 
Cy = Coo(1 + da*) (d constant for a given vertex angle), 
where 
C, = L/ipq?, W, Cp = D/3pq, W, Cy = M/[3pq W?, Coo = Cp for « = 0. 


The results obtained for 0C,,/0x, AC y,/0x, Cp) and d are tabulated in table 2. 
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These results are plotted in figures 11, 12 and 13. ‘The values of Cppo are 
compared with theoretical values given by Perry (1952) in figure 11. In 
order to compare the experimental values of the lift slope coefficient 0C;,/da, 











Figure 11. Measured drag coefficient as a function of o. 
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Figure 12. Measured lift slope coefficient as a function of o. 


with the theory developed in this paper, the results are extrapolated to give 
(0C,/ca) for ¢o=0. For values of B = 90°, 60° and 45°, there does not 
appear to be any dependence on cavitation number and accordingly the 
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results are averaged. For B = 30° and 15°, however, there is a significant 
increase of dC,/d« with the cavitation number oa, and a least squares fit 


is used to obtain (dC ,/0x), 


a 


The results for the moment coefficient are 
all fitted with a least squares fit to give (dC,,/d«), 


The values of 


V* “SN / . ° ~ 
(dC; /dx), _ and (dC3,/0x), _ , are compared with theory in figures 8 and 9. 
These results are tabulated in table 3. 
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Figure 13. Measured moment slope coefficient as a function of o. 
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Table 3. 


SSION OF RESULTS 


‘The drag results for zero incidence plotted in figure 11 show gocd 


ag 


agreement with theory. No theory is presently available for the dependence 


of the drag on the angle of incidence except for the case of a flat plate at 
zero cavitation number, when for a flat plate almost normal to the stream 


F.M. 


2S 
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the drag coefhicient may be expressed as C, = Cp, (1 — 0-000238x?) compared 
with an experimental value of Cp = Cp, (1 —0-0001712?) for x in the range 
O15 

Since the lift coefficient of a flat plate nearly normal to the stream may 
be expressed as Cptan« (neglecting skin friction), and as Cp, behaves as 
| +o, it would seem reasonable to suppose that the lift coefhicient would 
experimental results do not 


— “y 1 
Denave In a Similar manner. However, the 






show this trend, and the reason for this discre¢ yancy 1s not yet apparent. 
The extrapolated values of the lift slope coefficient have been plotted in 
theory when it will be seen that agreement is 
fair, part of the ditference between theory and experiment being accounted 
for by the fact that insufficient results were obtained to achieve great accuracy 
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APPENDIX I 
The evaluation of the integral I occurring in equation (13) 
Ixpanding in powers of a by ‘laylor’s theorem, and assuming that 
is of the same order of magnitude as a, we find that 
/ isin (Sem)[1 +e + Se? {h(1 — te) w(} —te)}]4 
+[—3aL,+pL,+3aL, +4 uL,| + O(a?), 
where (x) = (d/dx)log I(x), 


2 


L, cos acosea sina da, 


uv 


a+/9 


EB cos osin(1 —e)osin~’o do, 


L. e(1 — €*)z7/6, 
i, (1 —e€)?x/2. 
By integrating L, twice by parts, we find that 
L. cos(Se7) — (1 - e)*[1,_.—L,], 


where l, | sin n(sin #)-! dé, 


- cos($n7)[G(n) — G(—n)], 


with G(n) = 3[%(4 — jn) —b(3 — fn)]. 
Thus we see that 
| 1(1 —e)*x — cos(sem)[2 —€ +(1—€)?G(e)]. 
Similarly, integrating L, twice by parts, we find that 
L, - 2 cos(hem) + 4eL,— be(1 —¢)[/,—1,_,], 
so that 
L, = e(1 —&*)z/6— 4 cos(Sem)[2 — e? — e(1 — e?)G(e)]. 
Thus, on applying condition (6), we find that 
pe = at(e), 
-€?—¢(1—e?)G(e) 
e+(1 -€)2G(e) 


where 





This result is consistent with our original assumption that yu and a are 
of the same order of magnitude. 


APPENDIX I] 
(a) The calculation of h = h(n) 


We have, from equation (18), 





7{coshh— 1] = J,(n)—coshh-/,(n) —sinhhJ,(n), (11.1) 
where 
r*1+cost . ™1+ cost . 
J,(n) =| — sin nt dt, J ,(n) | —_— sinntcost dt, 
SI i 9 sini 
and *1+cost : 
aa Jn) = | — cos nt sint dt. 





fy ‘Sind 
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The substitution ¢ = 2 leads immediately to 
Sy(n) = 3[Lonza tonal, 
in the notation of Appendix I, and hence 
J(n)+J3(n) =Jy(n+1) = 3onss tonsil 


J(u) —J3(n) = J,(n— 1) = 3 [Lona + Lon-s]- 
Solving equation (II.1) for e”, we find the values given in table 4. 


(b) The dimensions of the subsidiary cavity 
We have, from equation (15), on writing 2%*2Sk!'*"Z= ze'P and 
7 = 7/2—-t, 
dZ 1+cost 
int : | 
— = ——— e' [1 —cos(t—th)]. 
dt sin t ( ) 
If we write Z = X +7Y, then X is the variable along the top face of the 
wedge and Y the variable perpendicular to it. 
‘, of Y is the dimension of the subsidiary 
cavity perpendicular to the top face of the wedge, and is given by 
m14 cost ‘ , ; 
y=} —— [sin nt (1 —costcoshh)—cosntsint sinh h] dt, 
J, wD 


Then the maximum value } 


where ¢,, is the value of t where the expression in square brackets vanishes. 
The variable XY has, for n > }, two turning points X, and X, and the 
dimension of the cavity along the top wedge face is |X,—X.|; while, for 
n ~ 1, there is one turning point 3, the dimension of the cavity then 
being X,. It is seen that Y,, X., and X, are given by 
ape. patets 1+cost eee 
Ais Aaj As | — [cosut (1—cost coshh) + sinnt sint sinhh] dt, 


Die sint 





the appropriate value of ¢ being taken as the upper limit. | 
Numerically we find the values in table 4, | 





0 | 1-00000 | 0) 2-000 
| 0-2 | 1-19825 | 0-343 | 1-925 
| 0-4 1-38439 | 0-693 | 1-767 | 
i} 06 | 1:54252 | 1-004 1-526 
| og | 1-65503 | 1-221 1-296 
| 1-0 | 1-70711 | 1-298 1-106 | 
| | 
} | 











Table 4. 


AppenpIx I{] 
lhe function p(e) occurring in equation (14) is defined by 


ple) = pile) + te)pole), 


\ here 


pile) = —5e/24—63/4—-€/6 + 3e?% d,(e)+€(1 —€*)a,(e) 
po(e) = 1/24+ €/6—3e/4+e8 2—€ a,(1—«)+(1—€)*a,(e), 








~] 
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wherein 
‘0 cos o| cos Ao— 1+ 3A? sino] 
( 


lo 





/ eo 
axA) = 3+ | send 
n/2 sin” o 


and 
cos o[sin Ao — A sino — $A{ 1 —A*}sin? o} j 
¢ 





d.(A) = A/3+ | sin’ o - 


Computed values of a,(e), a,(1—e) and d,(e) are given in table 5. 








| | 
| € | a.(€) | a,(1—e) d,(€) | 
PRs ee eo oe FR core 
| 0-0 | 0-3333 | 0-4299 0-0000 | 
| O-2 | 00-3386 | 0-4034 0-0544 | 
| 0-4 | 0-3537 | 0-3764 } 0-1122 | 
| 0-6 | 0:3764 | 0:3537 | 0-1766 | 
| 0-8 | 0-4034 0-3386 | 0-2498 | 
} 1:0 | 0-4299 0-3333 0-3333 | 
| | | | 





Table 5. 
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SUMMARY 
A comparison is made between various methods for the 
calculation of the pressure distribution on an elliptic cone in 
supersonic flow. In particular, a criterion first put forward by 
Van Dyke (1956) is reconsidered. ‘The paper is presented in 
two parts under separate authorship, taking the form of a discussion 
of some controversial aspects of this subject. 


Part I. (By R. V.-L.) 

Inarecent paper Van Dyke (1956) derived the second-order slender-body 
solution for an unyawed elliptic cone in supersonic flow. Considering his 
results to be exact, he used them as a criterion for comparing various 
approximations in compressible flow theory. Among these the method of 
linearized characteristics (Ferri 1951) was discussed, in particular, its 
application to the analysis of the flow field about non-circular cones. This 
method in general proceeds by linearizing the departure between the flow 
field of interest and some known non-uniform flow close to the actual one; 
along these lines Ferri treated the particular problem mentioned above 
by superposing perturbations on the axisymmetric solution, the disturbance 
velocity components being expanded in Fourier series in the polar angle @ 
up tothe tenth term. With this approach the entire flow field and the shock 
wave can be determined over a wide range of Mach numbers, in contrast 
to orthodox methods of linearization. 

Van Dyke based his appraisal of Ferri’s approximation on the two 
following numerical comparisons instituted for an elliptic cone of 3:1 axis 
ratio and area equivalent to a 10° circular cone, at M = v2: 

(1) The law expressing the slender-body pressure distribution was 
expanded in a Fourier series, and the error introduced by 
termination of this series at cos 106 was evaluated. 

(2) An approximation to the second-order slender-body solution was 
developed following Ferri’s procedure of expanding the velocity 
components in a Fourier series in 6, while linearizing with respect 
to deviation of the cross-section from circular, and calculating 
pressure coefficients from the full relation; the resulting pressures 
were compared with the exact values. 
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Van Dyke suggested that the large discrepancies between his complete and 
approximate results gave a measure of the errors to be encountered in 
Ferri’s method, supposedly due to termination of the series at cos 106 
and to linearization of the boundary conditions. Since good accuracy is 
to be expected from the second-order slender-body theory for the example 
in question, and since Ferri (1951) reported good agreement with 
experiments for a similar case, it seemed interesting to ascertain the origin 
of these discrepancies. ‘The matter has been the subject of private 
communications with Dr Van Dyke, who has provided many interesting 
comments and much information; these are gratefully acknowledged. 
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Figure 1. Comparison between various approximations to the pressure 
distribution around an elliptic cone. 


It is shown in the present note that satisfactory results can be obtained 
from the method of linearized characteristics, even when only linear terms 
are retained in the boundary conditions, provided that an area rule 
requirement is satisfied. The pressure distribution predicted by this 
procedure for the aforementioned elliptic cone at WM = v2 compares 
favourably with the second-order results (see figure 1); moreover, the 
experimental data presented in figure 3 of Van Dyke’s paper indicate that 
the actual values in the region of small 6 (tan@ = b/atany in the notation 
of his paper) should fall between the two predictions. ‘Thus both analyses 
exhibit approximately the same accuracy at the surface of the body, even 
though the particular example (a slender body at fairly low Mach number) 
represents an unfavourable test for Ferri’s method. It is also pertinent to 
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observe that Van Dyke’s solution is not uniformly valid throughout the 
flow field and fails at the shock, where conversely the accuracy of the 
linearized characteristics method is best since the perturbations there are 
smallest. 

The area-rule requirement for the expansion of the boundary conditions 
is justified by the following considerations. If, following Ferri, the shape 
of the body cross-section is expressed in polar coordinates (7, #, 0) as the 
Fourier series 

bk, = Wy 4 2 ws», cosné (1) 

n= 1 
with only linear terms retained in ¢,,,, and if %, refers to the basic circular 
cone, the velocity components at the surface of the actual body are given by 


‘ Urn 
De) F [ 2 > ib o( 2") cos nt |, 
isc i ! CnnJ vb 


a Ww ' 
w= 2(v.») > ntbs,,,(— sin 6. 
on=1 Unn PO) 


In the notation used by Ward (1955) the boundary conditions are 


dy gs, * dv 
“ ~ dr se dr = S'(z), (3) 


. 
v, a2 Jo 2, Jodz 





where v is the distance along the direction normal to the body contour 
in a plane perpendicular to the free-stream direction z, and v,, v_ are the 
velocity components in the v-direction and z-direction respectively. 
Within the approximation of the linearized characteristics method, we have 


: Urn 3 
Vy, = (T,o)y,, SIN Yop§ 14 = bol 2( 2) + cot tan Joos nO >, 
/ ] Cnn Vinh J | 
} = U ? 
Uz = (Troy, cos bo, 1+ > & | 2( 2) — tan boy feos nos, 
- Ob F | dic 5) | 
\ } l Onn y J 


Oh 


(4) 


Substitution of these expressions into (3) yields 
202 tan he, = S'(2). 


‘This result can alternatively be obtained by observing that in a neighbourhood 
O[z,,,] of the body surface the following quantities are small of the same 
order (the symbol O[ ] indicating order of magnitude, and s entropy) 


1 os 
O[z,,] = O[w] ol = = | Oly, »]. 


If terms O[%,] are neglected, the equation of motion in_ spherical 
coordinates reduces to the incompressible form 
ev ow 
div v = 2z7,+v, cot#e+ — + —~ = 0 
ous sin ys 0@ 


’ 


and the flow in the neighbourhood in question can be defined by an 
incompressible potential since the entropy is constant on the body and 
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the effect of ds/d% on the velocities can be neglected, being O[,?,]. 
Expansion of this potential in spherical harmonics and subsequent 
application of Gauss’s theorem to a control volume delimited by the actual 
body, by the basic axisymmetric body and by two spherical surfaces of 
arbitrary radius leads to the same requirement for the zeroth term as 
quoted above. ‘Therefore, the coefficients of the Fourier series representing 
the cross-section of the body within the approximation of the linearized 
characteristics method should be obtained by an expansion satisfying the 
requirement of equivalent areas. For the particular case of the elliptic 
cone at zero angle of attack, this can be obtained either (i) by determining b 
in the equation 
bab) 

% = tan! ——_,—__—______. . 
[(a/b)? sin? @ + cos? @]1” 
so as to satisfy the condition ztan?y,, = S (approximation A in figure 1), 
or (11) by determining % from the expansion 


? a*b ~ fa—b\" 
tan?ys = => ———. = ab} 14+2 Y ait cos 2n8 }, 
a® sin? 6 + b? cos? 4 = \a+b) 


which, within the linearized approximation, gives (approximation B in 
figure 1) 


: 0 fa—b\n 
tan? yy, = ab; ib = doy, +sindy, cospy, > () cos 2nd. 
neat Yosene 





‘The expansion used in method (11) has been suggested in a_ private 
communication by Van Dyke. 

The poor accuracy exhibited by the results of Van Dyke’s linear 
perturbation of a circular cone depended on the use of a direct Fourier 
expansion for the body geometry which yields %», = 9-27 instead of 
uo, = 10°; in the chosen conical reference system, higher-order terms 
in %,,, were thus considered for the boundary conditions, while the linearized 
approximation was retained for the computation of the velocity components. 
In this connection it is interesting to observe the close agreement over the 
range 20 6 < 90° between Van Dyke’s approximate results and the 
curve obtained by imposing his boundary conditions in the linearized 
characteristic analysis (curve C in figure 1). The increasing discrepancy 
between the aforesaid curves in the range 0 @ < 20° indicates that 
linear perturbations of the basic flow field lead to different accuracy in 
the two cases. ‘This behaviour may in turn be attributed to the fact that 
the basic flow used in the linearized characteristics approach represents 
a uniformly valid solution, in contrast to the basic flow considered in 
Van Dyke’s analysis. 

The present example of a linearized characteristics calculation has also 
been extended to include terms up to cos18@*. The results, which are 


* Tables of the disturbance velocities up to cos 18@ for a basic 10°-cone have 
been computed by the Research Department of the Grumman Aircraft Engineering 
Company. 
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included in figure 1, are practically identical to those obtained by the 
approximation as far as cos 106, and thus indicate that the question of 
convergence of the Fourier series in Ferri’s method ts related to the problem 
of representing the body cross-section rather than the problem of describing 
the perturbation flow field. It is believed that Van Dyke’s expansion of the 
slender-body pressure distribution to cos 106 afforded a lower accuracy than 
that obtained in the present analysis because the disturbance velocities 
considered in the two cases are of different order of magnitude; indeed, 
the latter expands the differences in velocity components (referred to 
conical coordinates) between the actual and the basic non-linear conical 
How, while the former involves the changes in velocity from free-stream 
magnitude and direction. 


Part II. (By M. D. V. D.) 


Professor Vaglio-Laurin makes the interesting suggestion that the range 
of applicability of Ferri’s method may be increased by modifying it in 
accord with the equivalence (or ‘area’) rule due to Oswatitsch (1952). 
The modification involves the same linearization as Ferri’s original method, 
differing from it only in higher-order terms; so that a demonstration of its 
superiority cannot be based on the linearized equations. Instead, one must 
compare with a more exact (non-linear) theory or with experiment, both 
of which will be done here. 

First, in keeping with the spirit of my original paper (Van Dyke 1956), 
the second-order slender-body solution derived there has been used as the 
basis for a test of Vaglio-Laurin’s proposal. ‘That is, the second-order 
slender-body approximation is taken as a model of the full solution, and the 
simplifications of Ferri’s method and Vaglio-Laurin’s modification are 
introduced as additional approximations. ‘The results are shown in figure 2 
(which is the left half of figure 5 of my previous paper with two added 
curves). Vaglio-Laurin’s approximation 4 is seen to give a definite 
improvement in the general level of pressure, though the peak is blurred 
(so that at this moderate Mach number simple slender-body theory is more 
accurate, as shown in figure 5 of my previous paper). 

Second, Jorgensen (1957) has recently measured the pressures over 
an elliptic cone close in shape to that just considered. It has the same axis 
ratio of 3:1; and although its cross-sectional area is that of a 7:76 rather 
than a 10° circular cone, it was tested at a Mach number of 1-97 rather 
than 1-41, so that the supersonic similarity parameter (the ratio of cone angle 
to free-stream Mach angle) is actually somewhat higher. ‘The results are 
shown in figure 3. Again the actual pressure peak is seen to be rounded off 
by approximation 4. 

Vaglio-Laurin’s last paragraph is quite correct. Curtailment of the 
Fourier series at cos 10@ introduces no appreciable error; and $3.5 of my 
previous paper was wrong on this point for just the reason that he suggests. 
Hence my previous assessment of Ferri’s method stands except that the 








Higher-order approximations for the flow fieid about a cone 643 





. 


f Model of exact solution 
(2nd-order slender-body approx.) 









Model of Vaglio-Laurin’s method: 
- approximation A 


a ee ae B 


(2) 











/ 
Model of Ferri's method— 














O 
: 45 90 
n, degrees 


Figure 2. Further comparison between approximations to the pressure 
distribution around an elliptic cone. 
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Figure 3, Comparison between theoretical and experimental results. 
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error is to be attributed almost entirely to non-linear effects rather than 
partly to the curtailment of the series. The method in its original form 
cannot reasonably be applied to cones whose cross-section differs from 
circular by as much as a 3:1 ellipse. For the modification suggested by 
Vaglio-Laurin, however, such eccentricity is seen to be perhaps tolerable 
at high supersonic Mach numbers where no more accurate theory exists. 
For less extreme eccentricities Vaglio-Laurin’s method will yield good 
accuracy. 

In view of this conclusion, the numerical results of Ferri, Ness & Kaplita 
(1953) must be regarded as unreliable except for nearly circular bodies. 
For example, the drag of an elliptic cone decreases much less with flattening 
than is suggested there. Indeed, Jorgensen’s measurements at Mach 
numbers of 1-97 and 2-94 show that the decrease in wave drag is so slight 
as to be entirely offset by the increase in skin friction due to greater 
surface area. 
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SUMMARY 
In a stationary weak sound wave the four gas dynamical 
quantities, entropy, stagnation enthalpy, mass flow per unit 
area and impulse per unit area are constant. ‘The six processes 
in which pairs of these four variables are kept constant are 
studied for the case of any single phase fluid or mixture of 
fluids in equilibrium, and it is shown that the remaining variables 


are stationary at sonic points. Such points are shown to occur 
once only in each process and are identified as maxima or minima, 
on the assumption that the fluid is a normal one which expands 
on heating at constant pressure and for which (07p/ov"), 1s 
positive. 

Newton’s theory of sound assumed that the fluid temperature 
was invariant. Across a stationary Newtonian sound wave, the 
four quantities, temperature, mass flow and impulse per unit 
area, and a dynamical variant of the Gibbs function are constant. 
‘he six processes in which pairs of these four variables are kept 
constant are studied, and it is shown that the remaining variables 
ire stationary at points where the fluid speed equals Newton’s 
sound velocity. ‘These points are shown to occur once only in 
each process and are identified as maxima or minima with the 
further proviso that (0*p/0v"),, is positive. 

The two sets of processes have one member in common, 
that usually referred to as the Rayleigh line. ‘The Fanno line 
and the usual isentropic ‘nozzle’ process also belong to the first 
SEt. 

Finally, the variation of stagnation temperature and pressure 
in some of the processes and in stationary shocks is investigated. 


1. INTRODUCTION 


Steady, one-dimensional gas dynamics has been analysed by many 
authors, particularly for the case of a perfect gas. A certain amount oft 
work has also been done for arbitrary, single-phase fluids in equilibrium, 
with properties limited only by stability and the laws of thermodynamics 
(Courant & Friedrichs 1948; Kline & Shapiro 1954). It has been shown, 
for instance, that the speed of sound is a critical velocity for isentropic 
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nozzle flows, and for Rayleigh and Fanno processes, whatever the fluid. 
This is too striking a fact to be pure coincidence. It is the aim of this paper 
to construct a symmetrical mathematical treatment of these three and 
related processes, showing the dominant role of the sound speed. 

The importance of studying arbitrary fluid dynamics is increasing 
since high-speed flows in which real gas behaviour, reaction, dissociation 
and ionization occur are becoming of practical interest. The theory of the 
arbitrary fluid in this and other papers is applicable to reactive gas mixtures 
provided only that at each section of flow investigated there is thermodynamic 
equilibrium. ‘The speed of sound occurring in the analysis refers to the 
speed of waves of low enough frequency for chemical equilibrium to prevail. 
‘The theory cannot be applied to detonation waves in pre-mixed gases since 
there is true equilibrium on only one side of these waves. ‘The flow behind 
and relative to a detonation wave can be sonic according to the Chapman-— 
Jouguet rule, whereas behind a shock wave in a reactive gas mixture, which 
is in equilibrium on both sides of the wave, the flow must be subsonic. 
This last fact is evident from the work of Kline & Shapiro (1954), provided 
the gas mixture’s properties are such that (0"p/dv?), and (dp/0s),, are constant 
in sign. p,wv ands are pressure, specific volume and entropy, respectively. 

Gases such as HCl which dissociate without total molal or volumetric 
change may be treated by perfect gas dynamical methods (Shapiro & 
Hawthorne 1947) with allowance for variation in effective specific heats. 


2. THE GAS DYNAMICAL RELATIONS 

Reactive or inert gas mixtures and single phase pure substances in 
equilibrium are usually characterized by the fact that any two independent 
properties define their thermodynamic state uniquely. As a generalization 
of this, the state of each section of flow in a one-dimensional gas dynamical 
process may be defined by three independent properties. One choice would 
be the velocity u together with any two independent thermodynamic ones. 
Another selection is Ff, G and H, in which 


F = p+pu*, the impulse per unit area, 
G = pu, the mass flow rate per unit area, and 
H =h+u?, the stagnation enthalpy. 


p and h denote density and specific enthalpy respectively. For convenience 
F is chosen equal to p+ pu? rather than the more usual quantity 


(p + pu”) x (duct area of cross-section). 


F/G represents this latter quantity per unit rate of mass flow. ‘The symbol .”” 
is chosen rather than /) to avoid complicated suffices. Note that F = p+Gu. 

Specified values of F’, G and H usually correspond to two possible states, 
such as occur on the two sides of a stationary shock wave. Across a 
stationary weak shock or acoustic wave, s is constant in addition to F, G 


and H, a very significant fact. 
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The familiar reversible, adiabatic ‘nozzle’ process and Rayleigh and 
Fanno processes each involve the invariance of two of the properties F, G, 
f1 and s. ‘Vhere are three other such processes. All six will be studied 
simultaneously. With suitable assertions regarding reversibility, adiabatic 
or other conditions, the six processes may be interpreted physically in 
various ways. Variation of G will be taken to correspond to duct area 
change, although it could equally well correspond to change in the mass 
flow due to injection or withdrawal of fluid. ‘Typical interpretations of 
the six processes appear in table 1. 


Constant | ‘Typical interpretation 
— — [ ee —— _ = - — | 
| Eos | Reversible, adiabatic ‘ nozzle’ process with area change. 
F,G | Rayleigh process. Frictionless, constant area process with heat | 
| exchange or Ohmic heating, etc. 
| G,H | Fanno process. Adiabatic, constant area process with force in 
| | direction of motion due to wall friction, wire mesh, electromagnetic 
| | forces in ionized gas, etc. No external energy exchange. 
G, s | Reversible, adiabatic work exchange in a duct of constant area by 
| actuator discs, electromagnetic induction, etc. 
Pes | Reversible, adiabatic process where work exchange and area change 
| combine to keep F' constant. 
| oh Adiabatic ‘ nozzle’ process in which F is constant owing to friction. 





Table 1. 


The last two cases are rather artificial, and are more interesting 
mathematically than practically. ‘here is difficulty in interpreting F 
when G is not constant. More artificial interpretations of the first four 
processes are also possible. For example, the Fanno process can represent 
a reversible process with suitably synchronized work and heat exchanges. 

There being only three independent properties of a flow-section, F, G, 
H and s must be related. ‘The differential equation expressing this relation 
is 

dk 4 pl'ds p dH+u dG, (1) 
which follows from the definition of the variables and from the equilibrium 
relation 7T'ds = dh—dp/p, T being the absolute temperature. Denoting 
the invariance of quantities by suffices, we then have 

dF,,=udGy,, etc. (2) 
(‘here are six such relations each involving two differentials. We shall 
proceed to relate all twelve differentials and show in particular that all 
vanish when the state is sonic. In a process with two of F, G, H and s 
constant, at points where the remaining two variables are stationary, the 
small change occurring between two closely adjacent flow-sections having 
the same values of F, G, H and s is virtually a stationary sound wave. 
Thus M, the local Mach number, must be unity. ™, or u/a, where a is 
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the isentropic sound speed , (¢p/dp),, is another gas dynamical property 
of each flow-section. A flow-section in a given state can be imagined as 
lying within any of innumerable processes, including the six currently 
under consideration. 

From (1), uv = (eF/eG),,, and p = (¢F/eH),,, and hence 


(0u/OH)g, = (0p/0G) ne. (3) 
here is a large array of these Maxwell-type relations, from which it is 


evident that all the differentials dH,,, dG,,,, etc., will vanish simultaneously 
| for flow-sections in certain states, provided there is no irregular behaviour 





of the quantities u, p, etc. Moreover, it is easily shown that 


dG ng p(1 . \[*) du, (4) 
so that sonic states are the relevant ones as expected. 

‘l'o provide a firmer basis for the analysis, it is convenient to take u as 
the major independent variable. Quantities of the type (CG/eu),,,, written 
G,,., will be studied and related. 

from (3) 
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It is necessary to study the behaviour 


nt thie quantities appearing } 
the table before the criticality of the sonic state can be clearly asserted. 
Vi? is well-behaved, as the speed 0 


ot sound must always be defined in real 
substances. ‘lhe quantity appearing in the (/, G) Rayleigh case 1s 


os c, [ / \ 
/ | = | = (10) 
( / a ( p } , 
which is well-behaved except at such rare points as water at + C at 
atmospheric pressure, where (¢ 


, ] 
(dl /cdp),. tends to 
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‘Table 2. Derivatives with resp ) ilti ( ] 
heat at constant pressure, must be positive stabiltt This is consist 
with the findings of Kline & Shapiro (1954) pointed out the critical 
of the sign of (cop os), and hence ot (¢s/cp ih yecific volum«e is 
reciprocal of p. ‘The quantity appearing in the (f°, 7)-case 1s 
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non-zero. ‘The quantity appearing in the (G,/7) bani is 
| ; 
joo \ ( / 19 
i | rs 
\ on} a*\ p)} \ 
is seen to be finite and positive except conceivably in a thermal contractio 
region, Dy reference to (11) above VV 


nay note that 





650 Ff. A. Sherchf 


which is the quantity studied by Kline & Shapiro (1954) in their discussion 
of Fanno processes. ‘They considered it very unlikely that 7 + 7(dp/ds), 
would ever be negative for real fluids. 

Summing up, we may say that if the above special states are excluded, 
the pair of the variables F, G, H and s that are not kept constant are 
stationary when and only when M = 1. 

‘The next question to consider is whether these stationary values ar 
maxima or minima. Kline & Shapiro (1954) have studied the (//,s), 
(F,G) and (G,H) cases and have shown separately in each case that the 
sign of (d%p/ev?). is the crucial factor. ‘lhis also decides whether com 
pressive acoustic waves steepen into shocks, and whether shocks must be 
compressive. 

The six cases, keeping pairs of F, G, H and s constant, will now be 
treated in parallel to determine the nature of the stationary values as w varies. 
lo fix the change of sign of F’, G’, H’ or s’ at the sonic state, it is necessary 
to study 0M?/du at the sonic state. It is found that ¢M?/cu takes the same 


value in all six cases when MW |. ‘This follows from the five equations: 


(2M? = ‘2 | 
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Equation (11) indicates that (¢p/cp), is not infinite, and so (da?/ées), is not 
ntinite 
We conclude that ¢M*/cu takes the same value when W 1, whichever 
pair of F, G, H and s is kept constant. ‘The common value is most easil\ 
evaluated in the (G,s) case, where w/v is constant 
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Kline & Shapiro (1954) indicated that (d?p/dv*), is positive for all single 
phase fluids for which accurate data are available. We therefore assum« 
that (0?p/dv*), is positive, and also that there is no thermal contraction, 
in what follows. It is then possible, using (10), (11) and (12), to deduce 
the nature of the stationary values of F, G, H and s at sonic points from 
the change in sign in the relevant derivative with respect to u. ‘Table 3 
lists the resuits. It should be observed that the sonic points have the same 


nature as for a perfect gas. 


Constant Nature of sonic points 
//, s (isentropic nozzle) F’ max.; G max. 
GJ, Hf (Fanno) F’ min.; s max. 
I’ G (Ravleigh) /1 max.; s max 
Gs F min.; H min 
Fy s (; max.; H min. 
1 Dd 5 | Gs max.; Ss max. 


Panle-3, 


If we assume that the variation of W, F, G, H and s with u in the above 
cases is continuous, it follows that only one sonic point occurs in each 
process (if (d"p/dv?), is fixed in sign) since maxima and minima must occur 
alternately. © (Sonic. points of inflexion would require 0M?/cu and 
(d2p/dv?), to vanish.) ‘This fact was deduced by Kline & Shapiro (1954) 
from a different argument in the Fanno and Rayleigh cases. 


3. ‘THE INFLUENCE OF THE SPEED OF NEWTONIAN SOUND 
Newton’s calculation of the speed of sound in air proved to be too low 
because of his assumption that sound propagation was an isothermal process 
instead of an isentropic one. Nevertheless the speed of Newtonian sound, 
which will be denoted by 7 to distinguish it from the adiabatic sound speed a, 
proves to have some physical significance. 
The symbol ' will be used for the ‘isothermal Mach number’, the 
ratio un. ‘Uhe specitic heat ratio y is equal to V?/M? because 
a® (dp/dop), (ds/eT), ae 
n (dplop), (0s/0T), Cy * 
y is always greater than unity since c,—c, = — T(évjeT){(cp/cv), and 
(6p/czv), must be negative for stability. y reaches unity only at points of 





zero thermal expansion where (ev eT), vanishes. In passing, it it worth 
noting that the speeds of adiabatic and Newtonian sound are (yp/p)'? and 
(p/p)'* for any fluid which obeys Boyle’s law, p/p = f(T). A reactive gas 
mixture whose mean molecular weight changes does not satisty this 
condition. 

It is well known that, for a perfect gas undergoing a Rayleigh process, 
the temperature has a stationary value when M = y"'* and N= 1. This, 

2T2 
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lowever, Is a general result, applicable to all single phase Huids it 


ur 3 and a? ~ ) ; 
\ ¢ p/y ( RP, 


It follows that (APS fi p 9 AT\ [ep \ 
ae a ’ M (=) (4) 


and hence that 
( 1 al 
M? = : 
l ( re ( os ) 


quilibrium, for 


(eT) 
\ as J 
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From this we have T 1-yM 
T= 
p 
vhich vanishes when M = y~-'*. If the fluid is a reactive gas mixture, 
must be taken as the effective specific heat ratio, allowing for reaction. 
\ flow-section such that NV = 1, i.e. capable of sustaining a weak steady 


Newtonian sound wave in which F, G and T do not change significantly, 
obviously corresponds to a state where 7 is stationary if F and G are fixed 
as in a Rayleigh process. It also corresponds to the section when G and 
the duct area have stationary values in a reversible, isothermal nozzle process. 
in fact there are again six simple processes for which n is the critical speed. 
\ unified, symmetrical treatment is again possible if a new property Z, 
a dynamic form of Gibbs function, is introduced such that 
Z=H-Ts = (+32, 

being the usual Gibbs function per unit mass, A— 7's. Z is not the 
stagnation Gibbs function, which is H—7)s, where 7), is the stagnation 
temperature. Z has no obvious physical significance except in tsothermal 
cases where its change measures reversible work exchanges (see $4). ‘Uhus 
Z is constant in the reversible, isothermal nozzle process. ‘Vhe lack of 
significance of Z is related to the fact that the datum for s ts arbitrary and 
therefore Z contains an arbitrary multiple of T when 7 varies. ‘The two 
cases where Z, but not 7, is constant will not be considered further. 


Constant ‘Typical interpretation 
F,G Ravleigh process, as already discussed 
‘Ree 2 Isothermal, reversible ‘ nozzle * process with area change. 
G, 1 Isothermal, constant area process with state change produced by wall 


friction or work exchange 

F, 7 Isothermal process in which area change and wall friction or work 
exchange combine to keep F' constant. 

lable 4. 

he four variables FY G, 7 and Z are related by the equation 


dF p dZ + ps dI'+u dG. (13 ) 


I'he four significant cases in which pairs of the variables are kept constant 


are interpreted in table +. 
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This situation may be analysed in a manner completely analogous to 
that in the foregoing section. For brevity, only the main results will be 
\gain uw is taken as the main independent variable and primes 


presented. 
‘Table 5 lists these derivatives as 


denote ditterentiation with respect to u 


multiples of (1 — V%). 


Constant G'/p F’ /pu 67’ fu Z fu 
Tee 
(isothermal | 9 v 
| ‘nozzle ’) 
| eal 
aTh 


(Rayleigh) 


7 (1+ N?) 0) 


‘Table 5. Derivatives with respect to v as multiples of (1 — N*) 


It is necessary to study the behaviour of the quantities appearin 
table 5 before the criticality of the Newtonian sonic state can be 
clearly asserted. N® is well behaved, since two-phase states in wl 
(Op/ep), and n® vanish are excluded. Also, 

(OT /0p),/ N? (OT /dp),)u? 


which is positive and finite if thermal contraction 1s excluded. Note that 


the factor s also occurs in the quantity s7”/u so that its sign ts irrelevant 
to the behaviour of 7’ in the Rayleigh case. 

The conclusion is that with the above restrictions the pair of the variables 
F, G, T and Z that are not kept constant are stationary when and only when 
V=1. ‘The nature of these stationary points depends upon the sign ol 


IN?/ou there. 
In the (G, 7’) case, it ts easily shown that 


IN2  v3/d2p) 
( 4 sihiiaii: Bie 14 
/ 


ov" 


Ou n? 2 
It is also easily demonstrated that 0N?/cu takes this same value at th 
The analvsis is similar to that 


Newtonian sonic point in all four cases. 
he condition 


in §2, with N,v, T, Z and ¢ replacing M, a, s, H and h. 
for 0N?/du to be identical is found to be that (en?/e7),, which is equal ti 
should not be 


v2(0*s/dv?),, and (Cn*?/cp),, which is equal to (07p/dp*),, 
infinite. No singularity upon isothermals such as would make the 
quantities infinite appears to be possible. We conclude that in all fou 
cases (14) holds. ‘The derivative (07p/év*), is usually positive, but cai 
thy 


become negative near the critical point for pure fluids. Tf we assum 
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the positive sign prevails, and that there is no thermal contraction, we can 
ascertain the nature of the stationary values when N= 1. ‘The results 
appear in table 6. ‘The behaviour is qualitatively like that of a perfect gas 
for substances for which there is no thermal contraction and for which 
(d%p/dv?),, is positive. 





: Nature of Newtonian sonic 
Constant 





points 
I’, G (Rayleigh) T max.; Z min. (if s 0) 
| T, Z (isothermal nozzle) F’ max.; G max. 
CR h F min.; Z min. 
jae de & G max.; Z max. 
| _ — J 
‘lable 6 
N= | 
oe 
oo SUBSONIC 
1M=| 
| * SUPERSONIC 
| / 
/ 








S 
Figure 1. he Rayleigh Line, for the case 


(ap év*), mr 0, (dp/0T), (), 


[f we assume that the variation of Ff, G, 7 and 7 with u in the above 
cases is continuous, it follows that only one point at which N = 1 occurs 
in each process (if (d*p/dv"), is fixed in sign) since maxima and minima 
must occur alternately and points of inflexion would require (07p/dv?),, to 


vanish. 

This, and the corresponding result for points at which .W/ = 1, con- 
strains a Rayleigh line on the (7, s)-plane always to have the form familiar 
for the perfect gas case, illustrated in figure 1. ‘I'he subsonic and supersonic 
branches cannot cross because this would imply two flow-sections with the 
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same values of 7’ and s, and hence in the same thermodynamic state with the 
same values of a and p, but with differing Mach numbers and velocities. 
But F = p+Gu, F and G being constant, which shows that for a given 
value of p, wu can have only one value. hus a crossing is precluded. 


4. SOME OTHER CONSIDERATIONS 
When propulsive effort is the prime concern, the quantity 
(p + pu”) x (duct area of cross-section) 

is of interest. ‘his quantity is represented by F/G per unit rate of mass 
How, 1/G being the area per unit rate of mass flow. It is possible to develop 
the analysis using the four variables F/G, 1/G, H and s, which are related 
by the equation 

ud( f'}G)+ T ds dH + pu d(1 G). 


‘The results differ little from those in § 2, there being only two fresh processes 
involved. ‘hese are the cases where F/G and H ors are constant. For these 
cases 1/G and s or H, respectively, are stationary at the sonic point. 

Another approach is offered by the selection of the four variables H, p, 
u and s tor analogous treatment. ‘lhe basic relation here is 


p dH = dp tpl'ds+G du. 


No critical point like the sonic point occurs in this case, however. for 
instance, when H and s are fixed, p is stationary only when uw and G vanish. 
Another case belonging to this family is that where s, p and hence the 
thermodynamic state are fixed. It corresponds to cases where the kinetic 
energy is converted reversibly into work by idealized impulse turbines or 
electromagnetic induction. Again the sonic point is not critical. 

There are obviously many other similar, alternative selections of variables 
possible. 

Interpretating the individual cases physically reveals an asymmetry in 
the way that heat and work influence the problem. In reversible cases, 
the external heat exchanged between two adjacent flow-sections is 7'ds 
per unit mass flow, but the corresponding work quantity cannot be expressed 
in terms of any single perfect differential in general. Alternative, general 
expressions for dW, the external work exchanged between two adjacent 
flow-sections per unit mass flow in a reversible process, are 


dW = Tds—dH = —(u du+dp/p) = (udG—dF)/p dZ—s aT. 


‘The difference between dW and the actual work exchanged in an irreversible 
case is the dissipation. In reversible cases without work exchange such 
as the two nozzle processes and the Rayleigh process, dW vanishes. ‘This 


will be recognized as indicating that Euler’s equation 
udu +-dp/p = 0 


is applicable. 








656 ¥. A. Sherchff 





Gas dynamical analysis frequently makes use of the stagnation state 
corresponding to a given flow-section. ‘his state, usually denoted by the 


suthx 0, is detined as that for which 
h, = H, Sy = S. 


he stagnation temperature 7), is of interest mainly for substances for which 


Uy 1), 
h is a tunction of temperature only, as for a perfect gas. For all normal 
Auids 7\, is larger than 7’, the ‘ static’ temperature, since H is always greater 
than A, and 


(< / l ( Op 
} 

ch | Yoana r) 

provided thermal contraction ts excluded. 


f 


dp, o,(a@H — T,, ds). 


Variation of the stagnation pressure p, 1s given by the relation 


p, is obviously constant in the isentropic nozzle process. In a Fanno 
process, H is constant and 


dp, = —py Ty ds. 


‘hus if the process is irreversible and adiabatic, p, must decrease. 
In a Rayleigh process, dH = 7'ds and 
dp Py \ [ 1’) as, { 15) 
hus the sonic point is the point of minimum p, tor normal fluids without 
thermal contraction and for which (0*p/dz?), is positive. 

In the process with s and F or G constant, py behaves qualitatively 
like H, reaching a minimum at the sonic point for normal fluids. In the 
(F, H) process it is again minimum when MV = 1 since s is then a maximum. 

‘hus, of the six F’, G, H, s cases, py is constant in one case and is minimal 
at the sonic point in all the others for normal fluids. 

In the isothermal nozzle process, dH = Tds and (15) again applies. 
Neither p, nor s has stationary points in this case, however. 

\cross a stationary shock, // is constant whereas s always increases. 
It follows that p alwavs decreases because (Op, Os) y : rived bie (), 


5. CONCLUDING REMARKS 


One conclusion of the analysis is to confirm and generalize the result 


t 
OT 


Kline & Shapiro (1954) that the gas-dynamical behaviour of single-phase 
fluids in equilibrium is qualitatively simular to that of a perfect gas, provided 
that there is no thermal contraction and that (0%p/dv). is positive. he 
idditional demand that (0?p/cv*), be positive is made. 


‘The speeds of adiabatic and Newtonian sound are seen to be significant 


ince a flow-section where the velocity is sonic is capable of undergoing 
an infinitesimal acoustic-type change in which F and G together with 
ither s and H (if W/ = 1) or 7 and Z (if N = 1) do not change. It is then 
ipperent, for instance, why 1 Ravyleich process with F and G constant, 


or 7 is stationary at the tw types Of SOIC pout 
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The tendency of nozzle, Rayleigh and Fanno processes to choke when 
the velocity becomes sonic is well known. ‘he similar choking occurring 
in the (G,s), reversible, adiabatic work-exchange process is perhaps a new 
result. ‘he minimum value of // which occurs at the sonic point means 
that a limit is placed on the work that may be extracted reversibly in a duct 
of constant cross-sectional area. ‘The comparable process mentioned in § 4, 
in which the duct area varies in such a way that work extraction changes 
only the kinetic energy and not the thermodynamical state of the fluid, shows 
no choking tendency. ‘lhe limit ou the work extraction here is provided 


by the growth in duct area as the velocity falls. 
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CORRIGENDUM 


“Wave tormation in laminar flow down an inclined plane’’, by 
‘T. Brooke BENJAMIN (7. Fluid Mech. 2, 1957, 554). 

Page 560. In equation (3.11) the equals sign following the first term 
should be replaced by a minus sign. 

Page 570. In equation (5.10) the numerical factor should be replaced 
by 0-448, and in equation (5.12) the numerical factor should be replaced 
by 0-00783. Equation (5.13) should read as 


J = exp(0-108R7%}. 
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REVIEWS 


Elements of Gasdynamics, by H. W. LikzpMANN and A. RosHko. 
New York: John Wiley & Sons, 1957. 439 pp. $11.00 or 88s. 


‘The term gas dynamics, which Liepmann and Roshko have chosen to 
combine into a single word, is now widely used to denote that area of fluid 
physics which deals with compressible flow and includes some of the more 
basic applications to aerodynamics. It is perhaps not more than a decade 
since the subject blossomed forth as a major category of aeronautical science, 
somewhat distinct from the more classical parts of fluid mechanics, and it is 
also a decade since the predecessor to the present book was published, namely 
Liepmann and Puckett’s Aerodynamics of a Compressible Fluid. In 1947 
and subsequent years, this work helped to fill an obvious need for an intro- 
ductory gas dynamics book written in English. It gave an admirable treat- 
ment of the fundamentals and showed a wise choice of applications. — Prior 
to 1947, coverage in the book literature was inadequate (with some 
exceptions such as Sauer’s book Theoretische Einftthrung in die Gasdynamik 
(1943)), which is curious considering that the subject has roots extending 
well into the past century ; for example, the work of St. Venant and Wantzel, 
Riemann, Mach and Rankine. Also, such investigations as Prandtl’s 
schlieren measurements of waves in supersonic channel flow (1907) and 
Meyer’s basic theoretical work on supersonic plane flows (1908) were four 
decades old at the end of World War II. 

There are two principal reasons why the writing of a gas dynamics book 
represents a greater challenge now than ten years ago. Firstly, a number of 
authoritative books on the subject have already appeared during the past 
decade, which give a rather wide coverage and diversity of emphasis among 
component topics. ‘The authors or editors, as the case may be, include 
Ferri, Courant and Friedrichs, Howarth, Oswatitsch, Shapiro, as well as 
Emmons, who edits the forthcoming gas dynamics volume of the Princeton 
series. Secondly, the discipline itself has been extended in scope. The 
subject matter has grown in volume and its conceptual boundaries have 
penetrated more deeply into the realm of particle physics, as energies and 
transit times have approached those characteristic of molecular processes. 
Thus, the problem of optimizing the choice of material for a book of reasonable 
size is more severe. According to the preface, it was the first of these reasons 
which led the authors to write a completely new book rather than revise the 
Liepmann- Puckett work. ‘The second reason generated the decision to 
divide the material into two independent volumes of which this is the first, 
with the second yet to be written. 

In this reviewer’s opinion, the authors have, in fact, implemented etfec- 
tively their aim of providing a working understanding of the essentials of 
high speed flow. ‘he book combines an approach characteristic of a first 


class textbook with a careful choice of subject matter based on considerable 
classroom experience. As such, it should be appreciated particularly by 
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new students of gas dynamics. On the other hand, much of the contents 
should be largely familiar to those who already claim some competence in the 
field, and such persons will find that the book is well suited as a reference and 
review text. 

The organization of material is largely parallel to that in Liepmann and 
Puckett. This is, if anything, another point in its favour. Part I of the 
older book has been replaced by the first six chapters of the present work and 
part [I by the last eight chapters. ‘The early chapters deal with one-dimen- 
sional flow where construction of the equations of momentum, continuity 
and energy is feasible in conjunction with simplified physical concepts. 
As is known, much ground can be covered before further generalization 
becomes necessary, and good advantage is taken of this opportunity here. A 
broader treatment of the basic equations takes place in two widely separated 
steps. ‘The first of these occurs around the middle of the book and deals with 
the equations of frictionless low. Some instructive features of the presenta- 
tion at this point are the introduction and use of tensor notation, and an 
explicit comparison between the control surface or integral method of deduc- 
ing the equations and the use of the Eulerian derivative. In connection with 
the tensor formulation, the reader may wonder whether the authors succumb 
to an aesthetic urge to use mathematical finery at the expense of physical 
concepts; he can be assured that this is not so, for Liepmann and Roshko 
are too much concerned with the clarity of the underlying physics. 

The second step toward generalization occurs near the end of the book 
and deals with the formulation of the complete Navier-Stokes equations. 
Here, the tensor notation introduced earlier is put to work, and by means of a 
skilful presentation, the authors rapidly produce these equations. ‘The 
treatment is somewhat more sophisticated here and, accordingly, the section 
in question is officially included among the more advanced ones which may be 
omitted at a first reading. ‘This is, perhaps, consistent with the authors’ 
inductive approach, at least in the sense that it progresses from one-dimen- 
sional gas dynamics to the introduction of the Navier-Stokes equations over a 
span of twelve chapters. It is possible, however, that the development of the 
basic physical theory would be strengthened by including a preliminary 
discussion of the Navier—Stokes equations at an earlier point. As it stands 
now, this material is somewhat removed from the main stream of the book. 
Even the bare outlines of a formulation from the basic stress—strain dynamics 
of a fluid element would provide an initial framework for subsequent 
special cases, such as the excursion into one-dimensional ideal gas flow. 
‘That this approach gives the reader a broad perspective of physical concepts 
at an early stage has already been shown by other authors such as Prandtl in 
his article “The Mechanics of Viscous Fluids’”’, Volume II1, Section G 
ot Aerodynamic Theory edited by Durand (Berlin: Springer, 1935). 

‘The many satisfying features of the book include two chapters dealing 
with the similarity rules of high-speed and transonic flow. ‘The former 
presents a unified discussion which treats regimes from subsonic through 
to hypersonic. It is limited to small perturbation theory but, as the 
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authors point out, this 1s not especially restrictive from a practical point of 
view since modern acrodynamics is concerned mainly with thin wings. In 
the chapter on transonic flow, good use is made of the potential for plausi- 
bility arguments inherent in this branch of the subject. ‘There is also a 
fruitful comparison between transonic flow over wedges and cones, 
a subject which has been of particular interest to the senior author in 
connection with his own experimental work in this area. 

‘Turning to the matter of high speed flows with viscosity and conductivity, 
we recall that Liepmann and Puckett present an abbreviated treatment, too 
much so in the light of recent extended interest in the subject of boundary 
layers and other varieties of dissipative flows. Boundary layer theory ts of 
engineering importance not only because it yields numerical results for 
friction drag and heat transfer but also, as the present authors emphasize, 
since it confirms that viscosity does not influence the pressure field in high 
Reynolds number flow past a slender body. In addition, the discipline 
harbours a variety of physical concepts whose explicit treatment should be of 
considerable pedagogical value, and nowhere is the interplay of dynamic 
behaviour and intrinsic physical properties in better evidence. ‘The 
reviewer has in mind, for example, such matters as the variation of shearing 
stress distribution with temperature dependence of viscosity, interpretation 
of Prandtl number, its role in boundary layer energetics, and the examination 
of boundary-layer type Hows as diffusion processes. 

The treatment of non-ideal gas flows in the present book ts, tn fact, another 
one of its strong points. ‘lhe boundary layer concept ts developed in a way 
calculated to give the reader a proper physical insight, and newer material 
presented includes recent measurements of turbulent boundary layers at 
high Mach numbers, as well as a discussion of Couette flow of a dissociating 
gas. 

The last mentioned topic reminds us that a growing variety of so-called 
‘real gas’ effects is now under investigation as a result of greater interest in 
high temperature gas dynamics. ‘Their proper incorporation as part of the 
formal analytical foundation of the subject is not at all complete, as yet, but 
the authors manage to supply some clear thinking in their discussion of the 
problem. It seemed to the reviewer that their remarks will be particularly 
instructive to those who have tended to equate real vs ideal gases with 
particle vs continuum gas theories. Dithculties of this type arise because 
a particle theory is used to define the nature of various physical properties of 
the gas, such as those concerned with the transport of momentum and heat. 
It is possible of course to account for these two in particular on the basis of 
bulk coefhcients which appear in the continuum equations. Also, the 
‘imperfections ’ of dissociation, ionization and condensation can be treated 
by acontinuum analysis. However, in the case of, say, viscous shear by itself, 
one can visualize a rather satisfactory macroscopic physical model, whereas 
proper physical insight as a prelude to correct continuum-type treatment of 
several other real gas phenomena does require an evaluation of the molecular 
behaviour of the gas. The intrinsic non-equilibrium character of the 
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effects in question does not make the analysis any easier. ‘Vhat the dynamics, 
and particularly the energetics of high temperature gas How, does rest on this 
more complex foundation ts a point which Liepmann and Roshko manage to 
transmit successfully to the reader. ‘lhe same cannot be said for other 
current books in gas dynamics. 

‘he authors stress the point, and correctly so, that the key to the relation- 
ship between the continuum formulation and that of the kinetic theory lies 
in the effect of fluctuations due to molecular structure on the ability of 
measuring instruments to obtain correct mean values. With this as the 
proper criterion, they suggest that it is possible to apply a continuum analysis 
even to rarefied gas flows in which the Knudsen number is appreciably 
larger than unity. ‘This is done by measuring for a sufficiently long period 
of time and applying the ergodic hypothesis, which establishes an equivalence 
between time and space averages. (‘The approach is not new to the turbu- 
lence investigator who performs a statistical evaluation of hot-wire fluctuation 
data.) Beyond a certain point, however, the molecular fluctuations would 
grow so large that it is not practical to attempt the correct measurement of 
mean values. ‘I'he kinetic theory treatment then becomes the proper one. 
The motivation on the part of some for stating unequivocally that the con- 
tinuum theory breaks down for rarefied gas flow probably stems, in the 
authors’ opinion, from attempting to use a solution of the Navier—Stokes 
equation, e.g. a high Reynolds number solution, which does not apply in the 
first place. 

Each succeeding book in gas dynamics these days appears to include a 
correspondingly larger chapter devoted to a discussion of methods of physical 
measurement. ‘This is entirely warranted as there can be little doubt that 
the developing nature of the science places ever-growing stress on the results 
of pertinent experiments. ‘lhe discussion here is meant to centre around 
experiments of the wind-tunnel tvpe. Except for omitting free flight 
techniques, this is not too restrictive, since shock tubes and other variations 
of the more conventional tunnel now require a wide variety of experimental 
techniques. ‘There is probably a need for a completely new book in this field 
to supplement the 1954 volume of the Princeton series, Physical Measurements 
tm Gas Dynamics and Combustion. An extended list of contents for such a 
work would, indeed, include several intriguing topics, such as recent experi- 
mental work in magneto-gasdynamics. ‘The admittedly limited treat- 
ment by Liepmann and Roshko is, nevertheless, a substantial one marked by 
an excellent choice of material which is well presented. It represents a 
most suitable introduction to experimental problems for the student. ‘This 
is hardly a surprise as both authors are first class experimentalists in their own 
right. 


I>. BERSHADER 
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